diff --git a/.gitignore b/.gitignore index 8dddad7..77ff7ad 100644 --- a/.gitignore +++ b/.gitignore @@ -4,7 +4,15 @@ scripts venv -build* +build/ +build-*/ compile_flags.txt ._* tensor_test_files + +# MATLAB/Octave MEX files +*.mex +*.mexa64 +*.mexw32 +*.mexw64 +*.mexmaci64 diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index b750526..6cb3e62 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -8,6 +8,7 @@ repos: rev: v16.0.6 hooks: - id: clang-format + exclude: '\.m$' - repo: https://github.com/pre-commit/pre-commit-hooks rev: v4.4.0 diff --git a/CMakeLists.txt b/CMakeLists.txt index 6ea7531..91247b2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,11 +9,18 @@ cmake_policy(SET CMP0079 NEW) set(CMAKE_C_STANDARD 11) set(CMAKE_CXX_STANDARD 20) -set(CMAKE_C_FLAGS "-O3 -march=native") + +# FIXME: -march=native is not portable +# set(CMAKE_C_FLAGS "-O3 -march=native") +set(CMAKE_C_FLAGS "-g ") +set(CMAKE_CXX_FLAGS "-g ") + +option(ENABLE_SANITIZERS "Enable Clang sanitizers" OFF) include(GNUInstallDirs) add_library(binsparse STATIC) +add_library(binsparse_dynamic SHARED) add_subdirectory(include) add_subdirectory(src) @@ -24,8 +31,15 @@ add_subdirectory(src) find_package(HDF5 REQUIRED COMPONENTS C) target_link_libraries(binsparse PUBLIC ${HDF5_C_LIBRARIES}) +target_link_libraries(binsparse_dynamic PUBLIC ${HDF5_C_LIBRARIES}) include(FetchContent) + +# Force cJSON to build as static library for MEX compatibility +set(BUILD_SHARED_LIBS_BACKUP ${BUILD_SHARED_LIBS}) +set(BUILD_SHARED_LIBS OFF CACHE INTERNAL "Build cJSON as static library") +set(ENABLE_CJSON_TEST OFF CACHE INTERNAL "Disable cJSON tests") + FetchContent_Declare( cJSON # GIT_REPOSITORY https://github.com/DaveGamble/cJSON.git @@ -34,8 +48,12 @@ FetchContent_Declare( ) FetchContent_MakeAvailable(cJSON) +# Restore original BUILD_SHARED_LIBS setting +set(BUILD_SHARED_LIBS ${BUILD_SHARED_LIBS_BACKUP}) + configure_file(${cJSON_SOURCE_DIR}/cJSON.h ${CMAKE_BINARY_DIR}/include/cJSON/cJSON.h) -target_link_libraries(${PROJECT_NAME} PUBLIC cjson) +target_link_libraries(binsparse PRIVATE cjson) +target_link_libraries(binsparse_dynamic PRIVATE cjson) # Set up include directories properly for both build and install target_include_directories(${PROJECT_NAME} @@ -46,6 +64,14 @@ target_include_directories(${PROJECT_NAME} $ ${HDF5_INCLUDE_DIRS}) +target_include_directories(binsparse_dynamic + PUBLIC + $ + $ + $ + $ + ${HDF5_INCLUDE_DIRS}) + # Installation rules - these are always needed when the library is built install(TARGETS binsparse EXPORT binsparse-targets @@ -53,6 +79,12 @@ install(TARGETS binsparse ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} INCLUDES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) +install(TARGETS binsparse_dynamic + EXPORT binsparse-targets_dynamic + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} + ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} + RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} + INCLUDES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) # Install headers install(DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/include/binsparse @@ -65,6 +97,10 @@ install(EXPORT binsparse-targets FILE binsparse-targets.cmake NAMESPACE binsparse:: DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/binsparse) +install(EXPORT binsparse-targets_dynamic + FILE binsparse-targets_dynamic.cmake + NAMESPACE binsparse:: + DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/binsparse) # Create and install package config files include(CMakePackageConfigHelpers) @@ -84,6 +120,14 @@ install(FILES DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/binsparse) if(CMAKE_PROJECT_NAME STREQUAL PROJECT_NAME) + if (ENABLE_SANITIZERS) + set(SANITIZER_FLAGS "-fsanitize=address,undefined") + target_compile_options(binsparse INTERFACE ${SANITIZER_FLAGS} -g -O1 -fno-omit-frame-pointer) + target_link_options(binsparse INTERFACE ${SANITIZER_FLAGS}) + target_compile_options(binsparse_dynamic INTERFACE ${SANITIZER_FLAGS} -g -O1 -fno-omit-frame-pointer) + target_link_options(binsparse_dynamic INTERFACE ${SANITIZER_FLAGS}) + endif() + add_subdirectory(examples) add_subdirectory(test) endif() diff --git a/README.md b/README.md index a806025..29a5a82 100644 --- a/README.md +++ b/README.md @@ -18,7 +18,8 @@ This library provides a C interface for reading and writing binsparse matrices. #include int main(int argc, char** argv) { - bsp_matrix_t mat = bsp_read_matrix("chesapeake.bsp.hdf5"); + bsp_matrix_t mat; + bsp_read_matrix(&mat, "chesapeake.bsp.hdf5", NULL); if (mat.format == BSP_COO) { float* values = mat.values.data; @@ -33,6 +34,7 @@ int main(int argc, char** argv) { bsp_get_matrix_format_string(mat.format)); } + bsp_destroy_matrix_t(&mat); return 0; } ``` diff --git a/bindings/matlab/Contents.m b/bindings/matlab/Contents.m new file mode 100644 index 0000000..015bf6b --- /dev/null +++ b/bindings/matlab/Contents.m @@ -0,0 +1,35 @@ +% Binsparse MATLAB/Octave bindings +% +% Core functions (MEX, with .m help stubs): +% binsparse_read - read a matrix from a Binsparse HDF5 file +% binsparse_write - write a matrix to a Binsparse HDF5 file +% binsparse_from_ssmc - convert SSMC A+Zeros to a Binsparse matrix struct +% binsparse_minimize_types - minimize value/index types in a Binsparse struct +% binsparse_write_string_dataset - write an HDF5 UTF-8 string dataset +% write_binsparse_from_matlab - write an SSMC Problem struct (delegates to generate_bsp_from_ssmc) +% +% MATLAB helpers: +% generate_bsp_from_ssmc - write an SSMC Problem to a Binsparse file +% convert_to_problem_struct - convert a Binsparse problem to an SSMC Problem +% bsp_matrix_create - create a Binsparse matrix struct +% bsp_matrix_info - display information about a Binsparse matrix struct +% +% Build scripts: +% build_matlab_bindings - build all MEX functions with MATLAB's mex +% build_octave_bindings - build all MEX functions with Octave's mkoctfile +% compile_octave.sh - build the Octave MEX functions from the shell +% +% Tests: +% test_binsparse_read - error-handling tests for binsparse_read +% test_binsparse_write - write and round-trip tests for binsparse_write +% test_binsparse_from_ssmc - basic test for binsparse_from_ssmc +% test_binsparse_minimize_roundtrip - SSMC conversion + type minimization test +% test_bsp_matrix_struct - tests for the Binsparse matrix struct helpers +% test_convert_to_problem_struct - tests for Binsparse Problem conversion +% test_generate_bsp_from_ssmc - end-to-end test for generate_bsp_from_ssmc +% test_write_binsparse_from_matlab - end-to-end test for the SSMC writer MEX +% test_binsparse_roundtrip_dir - round-trip every .h5 file in a directory + +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause diff --git a/bindings/matlab/README.md b/bindings/matlab/README.md new file mode 100644 index 0000000..e09a708 --- /dev/null +++ b/bindings/matlab/README.md @@ -0,0 +1,315 @@ + + +# Binsparse MATLAB/Octave Bindings + +This directory contains MATLAB and GNU Octave MEX bindings for the Binsparse C library, providing a simple interface to access Binsparse functionality from both MATLAB and Octave. + +## Quick Start + +### Prerequisites + +**For MATLAB:** +1. **MATLAB** with MEX compiler support +2. **Binsparse C library** headers (included in this repository) + +**For Octave:** +1. **GNU Octave** with mkoctfile (usually included) +2. **C compiler** (gcc recommended) +3. **Binsparse C library** headers (included in this repository) + +### Setup MEX Compiler (if needed) + +**For MATLAB:** +If you haven't configured a MEX compiler yet: + +```matlab +mex -setup +``` + +Choose a compatible C compiler when prompted. + +**For Octave:** +Octave usually comes with mkoctfile pre-configured. Verify it works: + +```bash +mkoctfile --version +``` + +### Build the Bindings + +#### Option 1: MATLAB + +1. Open MATLAB and navigate to this directory: + ```matlab + cd('path/to/binsparse-reference-c/bindings/matlab') + ``` + +2. Build the MEX functions: + ```matlab + build_matlab_bindings() + ``` + +3. Test the installation: + ```matlab + test_binsparse_read() + test_binsparse_write() + test_write_binsparse_from_matlab() + ``` + +#### Option 2: Octave (from within Octave) + +1. Open Octave and navigate to this directory: + ```octave + cd('path/to/binsparse-reference-c/bindings/matlab') + ``` + +2. Build the MEX functions: + ```octave + build_octave_bindings() + ``` + +3. Test the installation: + ```octave + test_binsparse_read() + test_binsparse_write() + test_write_binsparse_from_matlab() + ``` + +#### Option 3: Octave (from command line) + +1. Navigate to this directory: + ```bash + cd path/to/binsparse-reference-c/bindings/matlab + ``` + +2. Build using the shell script: + ```bash + ./compile_octave.sh + ``` + +3. Test the installation: + ```bash + octave --eval "test_binsparse_read()" + octave --eval "test_binsparse_write()" + octave --eval "test_write_binsparse_from_matlab()" + ``` + +## Usage Examples + +### Reading Binsparse Files + +**In MATLAB or Octave:** + +```matlab +% Read a Binsparse matrix file +matrix = binsparse_read('path/to/matrix.bsp.h5'); + +% Read from a specific group +matrix = binsparse_read('path/to/matrix.bsp.h5', 'group_name'); + +% Matrix will be a struct with fields: +% - values: array of matrix values +% - indices_0, indices_1: row/column indices +% - pointers_to_1: pointer array (for CSR/CSC formats) +% - nrows, ncols, nnz: matrix dimensions +% - is_iso: boolean for iso-value matrices +% - format: string ('COO', 'CSR', 'CSC', etc.) +% - structure: string ('general', 'symmetric', etc.) +``` + +### Writing Binsparse Files + +```matlab +% Create a matrix struct (example: 3x3 COO matrix) +matrix = struct(); +matrix.values = [1.0; 2.0; 3.0]; +matrix.indices_0 = uint64([0; 1; 2]); % 0-based row indices +matrix.indices_1 = uint64([0; 1; 2]); % 0-based col indices +matrix.pointers_to_1 = uint64([]); % Empty for COO +matrix.nrows = 3; +matrix.ncols = 3; +matrix.nnz = 3; +matrix.is_iso = false; +matrix.format = 'COO'; +matrix.structure = 'general'; + +% Write to file +binsparse_write('output.bsp.h5', matrix); + +% Write with optional parameters +binsparse_write('output.bsp.h5', matrix, 'my_group'); +binsparse_write('output.bsp.h5', matrix, 'my_group', '{"author": "me"}'); +binsparse_write('output.bsp.h5', matrix, 'my_group', '{"author": "me"}', 6); +``` + +### Writing from SuiteSparse Matrix Collection Format + +```matlab +% Create a SuiteSparse Matrix Collection problem struct +Problem = struct(); +Problem.name = 'test_matrix'; +Problem.A = sparse([1 2 3], [1 2 3], [1.5 2.5 3.5], 4, 4); +Problem.title = 'Test Matrix'; +Problem.kind = 'artificial/test'; + +% Write directly from SuiteSparse format to Binsparse +write_binsparse_from_matlab(Problem, 'output.bsp.h5'); + +% Write with an explicit sparse format ('COO', 'COOR', 'CSC', or 'CSR') +write_binsparse_from_matlab(Problem, 'output.bsp.h5', 'CSC'); + +% Write with a format and gzip compression level (0-9) +write_binsparse_from_matlab(Problem, 'output.bsp.h5', 'COO', [], 6); +``` + +### Error Handling + +The MEX functions include proper error handling: + +```matlab +try + matrix = binsparse_read('nonexistent_file.bsp.h5') +catch ME + fprintf('Error: %s\n', ME.message) +end +``` + +## Files Description + +| File | Description | +|------|-------------| +| `binsparse_read.c` | MEX function for reading Binsparse matrix files | +| `binsparse_write.c` | MEX function for writing Binsparse matrix files | +| `binsparse_from_ssmc.c` | MEX function converting SuiteSparse A+Zeros to a Binsparse struct | +| `binsparse_minimize_types.c` | MEX function minimizing value/index types in a Binsparse struct | +| `binsparse_write_string_dataset.c` | MEX function writing HDF5 UTF-8 string datasets | +| `write_binsparse_from_matlab.c` | MEX entry point for writing an SSMC Problem (delegates to `generate_bsp_from_ssmc.m`) | +| `matlab_bsp_helpers.h` | Shared MATLAB/Binsparse conversion helpers for the MEX sources | +| `generate_bsp_from_ssmc.m` | Write a full SSMC Problem struct to one Binsparse file | +| `convert_to_problem_struct.m` | Convert Binsparse data back to an SSMC Problem struct | +| `bsp_matrix_create.m` | Utility function for creating matrix structs | +| `bsp_matrix_info.m` | Utility function for displaying matrix information | +| `build_matlab_bindings.m` | Build script for all MATLAB MEX functions | +| `build_octave_bindings.m` | Build script for all Octave MEX functions | +| `compile_octave.sh` | Shell script for building Octave MEX functions | +| `test_binsparse_read.m` | Test script for read functionality | +| `test_binsparse_write.m` | Test script for write functionality | +| `test_binsparse_from_ssmc.m` | Test script for SSMC conversion | +| `test_binsparse_minimize_roundtrip.m` | Test script for type minimization | +| `test_bsp_matrix_struct.m` | Test script for the matrix struct helpers | +| `test_convert_to_problem_struct.m` | Test script for Problem conversion | +| `test_generate_bsp_from_ssmc.m` | End-to-end test for the SSMC writer | +| `test_write_binsparse_from_matlab.m` | Test script for the SSMC writer MEX | +| `test_binsparse_roundtrip_dir.m` | Round-trip every .h5 file in a directory | +| `Contents.m` | Directory listing for MATLAB's `help` | +| `README.md` | This documentation file | + +## Technical Details + +### MEX Function Structure + +The Binsparse MEX functions demonstrate: + +1. **Header inclusion**: Proper inclusion of `` +2. **Type conversion**: Complete mapping between MATLAB and Binsparse types +3. **Error handling**: Using Binsparse error types (`bsp_error_t`) +4. **Memory management**: Safe allocation and cleanup +5. **MATLAB interface**: Proper MEX function structure with validation + +### Build Process + +**MATLAB build process:** + +1. Locates Binsparse include directory (`../../include/`) +2. Compiles MEX functions using MATLAB's `mex` command +3. Links against MATLAB MEX libraries +4. Validates compilation with test functions + +**Octave build process:** + +1. Locates Binsparse include directory (`../../include/`) +2. Compiles MEX functions using `mkoctfile --mex` +3. Links against Octave MEX libraries +4. Validates compilation with test functions + +## Extending the Bindings + +To add new Binsparse functionality: + +1. Create a new `.c` file with MEX function structure +2. Include `` and relevant headers +3. Add the filename to `mex_files` list in `build_matlab_bindings.m` +4. Create corresponding test functions + +### Example MEX Function Template + +```c +#include "mex.h" +#include + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { + // Input validation + if (nrhs != 1) { + mexErrMsgIdAndTxt("BinSparse:InvalidInput", "Expected 1 input argument"); + } + + // Call Binsparse functions + bsp_error_t error = /* your_binsparse_function() */; + + // Handle errors + if (error != BSP_SUCCESS) { + mexErrMsgIdAndTxt("BinSparse:Error", "%s", bsp_get_error_string(error)); + } + + // Return results to MATLAB + plhs[0] = /* create_matlab_output() */; +} +``` + +## Troubleshooting + +### Common Issues + +1. **MEX compiler not found** + - **MATLAB**: Run `mex -setup` to configure a compiler + - **Octave**: Install mkoctfile (usually comes with Octave) + - Ensure you have a compatible C compiler installed + +2. **Include paths not found** + - Verify you're running from the `bindings/matlab` directory + - Check that `../../include/binsparse/binsparse.h` exists + +3. **Compilation errors** + - **MATLAB**: Try building with verbose output: `build_matlab_bindings('verbose')` + - **Octave**: Try building with verbose output: `build_octave_bindings('verbose')` or `./compile_octave.sh --verbose` + - Check compiler compatibility with your MATLAB/Octave version + +### Platform-Specific Notes + +- **Windows**: + - MATLAB: Microsoft Visual Studio or compatible compiler + - Octave: MinGW-w64 (often included with Octave installer) +- **macOS**: Xcode command line tools required for both MATLAB and Octave +- **Linux**: GCC or compatible compiler should work for both MATLAB and Octave + +## Development Status + +This provides complete MATLAB/Octave bindings for Binsparse. Currently implemented: + +- ✅ Matrix reading (`binsparse_read`) +- ✅ Matrix writing (`binsparse_write`) +- ✅ Complete type support (all Binsparse types including complex numbers) +- ✅ Optional parameters (groups, JSON metadata, compression) +- ✅ Comprehensive error handling +- ✅ Build system and testing framework +- ✅ Round-trip compatibility (read → write → read) +- ⏳ Tensor support (future work) +- ⏳ Advanced sparse matrix operations (future work) + +## License + +This code is licensed under the BSD-3-Clause license, same as the main Binsparse project. diff --git a/bindings/matlab/binsparse_from_ssmc.c b/bindings/matlab/binsparse_from_ssmc.c new file mode 100644 index 0000000..c91722e --- /dev/null +++ b/bindings/matlab/binsparse_from_ssmc.c @@ -0,0 +1,388 @@ +/* + * SPDX-FileCopyrightText: 2024 Binsparse Developers + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +/** + * binsparse_from_ssmc.c - Convert SuiteSparse matrix A + Zeros to a Binsparse + * MATLAB struct. + * + * Usage in MATLAB/Octave: + * matrix = binsparse_from_ssmc(A) + * matrix = binsparse_from_ssmc(A, Zeros) + * matrix = binsparse_from_ssmc(A, Zeros, format) + * matrix = binsparse_from_ssmc(A, format) + * + * Inputs: + * A - sparse or dense matrix (MATLAB) + * Zeros - optional sparse matrix representing explicit zeros (same size as + * A) format - optional string: default 'COO' for sparse, 'DMAT'/'DVEC' for + * dense + * + * Output: + * MATLAB struct compatible with binsparse_write. + */ + +#include "mex.h" +#include +#include +#include + +#include "matlab_bsp_helpers.h" + +static bool array_uses_allocator(bsp_array_t array, bsp_allocator_t allocator) { + if (array.size == 0 || array.data == NULL) { + return true; + } + return array.allocator.malloc == allocator.malloc && + array.allocator.free == allocator.free; +} + +static bool matrix_uses_allocator(const bsp_matrix_t* matrix, + bsp_allocator_t allocator) { + return array_uses_allocator(matrix->values, allocator) && + array_uses_allocator(matrix->indices_0, allocator) && + array_uses_allocator(matrix->indices_1, allocator) && + array_uses_allocator(matrix->pointers_to_1, allocator); +} + +static bsp_error_t construct_array_with_allocator(bsp_array_t* array, + size_t size, bsp_type_t type, + bsp_allocator_t allocator) { + if (size == 0) { + array->data = NULL; + array->size = 0; + array->type = type; + array->allocator = allocator; + return BSP_SUCCESS; + } + return bsp_construct_array_t_allocator(array, size, type, allocator); +} + +static bsp_type_t sparse_value_type(const matlab_csc_t* matrix) { + return matrix->is_complex ? BSP_COMPLEX_FLOAT64 : BSP_FLOAT64; +} + +static void write_sparse_value(const matlab_csc_t* matrix, mwIndex src, + bsp_matrix_t* out, uint64_t dst) { + if (matrix->is_complex) { + double _Complex* out_values = (double _Complex*) out->values.data; + double imag = matrix->imag_values ? matrix->imag_values[src] : 0.0; + out_values[dst] = matrix->values[src] + imag * I; + } else { + double* out_values = (double*) out->values.data; + out_values[dst] = matrix->values[src]; + } +} + +static void write_explicit_zero(bsp_matrix_t* out, uint64_t dst) { + if (out->values.type == BSP_COMPLEX_FLOAT64) { + double _Complex* out_values = (double _Complex*) out->values.data; + out_values[dst] = 0.0 + 0.0 * I; + } else { + double* out_values = (double*) out->values.data; + out_values[dst] = 0.0; + } +} + +static void build_csc_merged(const matlab_csc_t* a, const matlab_csc_t* z, + bsp_matrix_t* out) { + bsp_error_t error; + + bsp_construct_default_matrix_t_allocator(out, bsp_matlab_allocator); + out->nrows = a->nrows; + out->ncols = a->ncols; + out->nnz = a->nnz + z->nnz; + out->format = BSP_CSC; + out->structure = BSP_GENERAL; + out->is_iso = false; + + error = construct_array_with_allocator( + &out->values, out->nnz, sparse_value_type(a), bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to allocate values array"); + } + + error = construct_array_with_allocator(&out->indices_1, out->nnz, BSP_UINT64, + bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to allocate indices array"); + } + + error = construct_array_with_allocator(&out->pointers_to_1, out->ncols + 1, + BSP_UINT64, bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to allocate pointers array"); + } + + uint64_t* out_colptr = (uint64_t*) out->pointers_to_1.data; + uint64_t* out_rowind = (uint64_t*) out->indices_1.data; + + out_colptr[0] = 0; + for (mwIndex j = 0; j < a->ncols; j++) { + mwIndex a_count = a->colptr[j + 1] - a->colptr[j]; + mwIndex z_count = z->colptr[j + 1] - z->colptr[j]; + out_colptr[j + 1] = out_colptr[j] + a_count + z_count; + } + + for (mwIndex j = 0; j < a->ncols; j++) { + mwIndex a_ptr = a->colptr[j]; + mwIndex a_end = a->colptr[j + 1]; + mwIndex z_ptr = z->colptr[j]; + mwIndex z_end = z->colptr[j + 1]; + uint64_t out_ptr = out_colptr[j]; + + while (a_ptr < a_end || z_ptr < z_end) { + if (z_ptr >= z_end || + (a_ptr < a_end && a->rowind[a_ptr] < z->rowind[z_ptr])) { + write_sparse_value(a, a_ptr, out, out_ptr); + out_rowind[out_ptr] = (uint64_t) a->rowind[a_ptr]; + a_ptr++; + } else if (a_ptr >= a_end || + (z_ptr < z_end && z->rowind[z_ptr] < a->rowind[a_ptr])) { + write_explicit_zero(out, out_ptr); + out_rowind[out_ptr] = (uint64_t) z->rowind[z_ptr]; + z_ptr++; + } else { + mexErrMsgIdAndTxt("BinSparse:DuplicateIndex", + "Duplicate indices between A and Zeros"); + } + out_ptr++; + } + + if (out_ptr != out_colptr[j + 1]) { + mexErrMsgIdAndTxt("BinSparse:InternalError", + "Merged column counts do not match"); + } + } +} + +static bsp_matrix_format_t parse_format(int nrhs, const mxArray* prhs[]) { + const mxArray* format_arg = NULL; + if (nrhs >= 3 && !mxIsEmpty(prhs[2])) { + format_arg = prhs[2]; + } else if (nrhs == 2 && mxIsChar(prhs[1])) { + format_arg = prhs[1]; + } + + if (!format_arg) { + return BSP_COO; + } + + if (!mxIsChar(format_arg)) { + mexErrMsgIdAndTxt("BinSparse:InvalidFormat", "Format must be a string"); + } + + char* format_str = mxArrayToString(format_arg); + if (!format_str) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", "Failed to read format string"); + } + + bsp_matrix_format_t format = bsp_get_matrix_format(format_str); + mxFree(format_str); + + if (format != BSP_CSC && format != BSP_CSR && format != BSP_COO && + format != BSP_COOR && format != BSP_DMAT && format != BSP_DMATC && + format != BSP_DVEC) { + mexErrMsgIdAndTxt("BinSparse:InvalidFormat", + "Supported formats: CSC, CSR, COO, DMAT, DMATC, DVEC"); + } + + return format; +} + +static void build_csc_from_a(const matlab_csc_t* a, bsp_matrix_t* out) { + bsp_error_t error; + + bsp_construct_default_matrix_t_allocator(out, bsp_matlab_allocator); + out->nrows = a->nrows; + out->ncols = a->ncols; + out->nnz = a->nnz; + out->format = BSP_CSC; + out->structure = BSP_GENERAL; + out->is_iso = false; + + error = construct_array_with_allocator( + &out->values, out->nnz, sparse_value_type(a), bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to allocate values array"); + } + + error = construct_array_with_allocator(&out->indices_1, out->nnz, BSP_UINT64, + bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to allocate indices array"); + } + + error = construct_array_with_allocator(&out->pointers_to_1, out->ncols + 1, + BSP_UINT64, bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to allocate pointers array"); + } + + uint64_t* out_colptr = (uint64_t*) out->pointers_to_1.data; + uint64_t* out_rowind = (uint64_t*) out->indices_1.data; + + for (size_t i = 0; i < out->nnz; i++) { + write_sparse_value(a, i, out, (uint64_t) i); + out_rowind[i] = (uint64_t) a->rowind[i]; + } + + for (size_t i = 0; i < out->ncols + 1; i++) { + out_colptr[i] = (uint64_t) a->colptr[i]; + } +} + +static void build_dense_matrix(const mxArray* mx_a, bsp_matrix_t* out, + bsp_matrix_format_t format) { + if (!mxIsNumeric(mx_a)) { + mexErrMsgIdAndTxt("BinSparse:InvalidMatrix", "Dense input must be numeric"); + } + + bool is_vector = (mxGetM(mx_a) == 1) || (mxGetN(mx_a) == 1); + + if (is_vector && format != BSP_DVEC) { + mexErrMsgIdAndTxt("BinSparse:InvalidFormat", + "Dense vector requires DVEC format"); + } + + if (!is_vector && format != BSP_DMAT && format != BSP_DMATC) { + mexErrMsgIdAndTxt("BinSparse:InvalidFormat", + "Dense matrix requires DMAT or DMATC format"); + } + + bsp_construct_default_matrix_t_allocator(out, bsp_matlab_allocator); + out->format = format; + out->structure = BSP_GENERAL; + out->is_iso = false; + + size_t nrows = mxGetM(mx_a); + size_t ncols = mxGetN(mx_a); + size_t total = mxGetNumberOfElements(mx_a); + + if (is_vector) { + out->nrows = total; + out->ncols = 1; + } else { + out->nrows = nrows; + out->ncols = ncols; + } + + out->nnz = total; + + bsp_error_t error = + matlab_to_bsp_array_allocator(mx_a, &out->values, bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to allocate dense values"); + } +} + +void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { + if (nrhs < 1 || nrhs > 3) { + mexErrMsgIdAndTxt( + "BinSparse:InvalidArgs", + "Usage: matrix = binsparse_from_ssmc(A [, Zeros] [, format])"); + } + + if (nlhs > 1) { + mexErrMsgIdAndTxt("BinSparse:TooManyOutputs", "Too many output arguments"); + } + + const mxArray* mx_a = prhs[0]; + bsp_matrix_format_t target_format = parse_format(nrhs, prhs); + + if (!mxIsSparse(mx_a)) { + bsp_matrix_t dense_matrix; + build_dense_matrix(mx_a, &dense_matrix, target_format); + plhs[0] = bsp_matrix_to_matlab_struct(&dense_matrix); + bsp_destroy_matrix_t(&dense_matrix); + return; + } + + if (target_format == BSP_DMAT || target_format == BSP_DVEC) { + mexErrMsgIdAndTxt("BinSparse:InvalidFormat", + "Sparse matrix cannot use DMAT/DVEC formats"); + } + + if (!mxIsDouble(mx_a)) { + mexErrMsgIdAndTxt("BinSparse:InvalidMatrix", + "A must be a sparse double matrix"); + } + + const mxArray* mx_zeros = NULL; + if (nrhs >= 2 && mxIsSparse(prhs[1])) { + mx_zeros = prhs[1]; + } + + if (mx_zeros && (mxIsComplex(mx_zeros) || + (!mxIsDouble(mx_zeros) && !mxIsLogical(mx_zeros)))) { + mexErrMsgIdAndTxt("BinSparse:InvalidZeros", + "Zeros must be a real sparse double or logical matrix"); + } + + if (mx_zeros && + (mxGetM(mx_a) != mxGetM(mx_zeros) || mxGetN(mx_a) != mxGetN(mx_zeros))) { + mexErrMsgIdAndTxt("BinSparse:DimensionMismatch", + "A and Zeros must have matching dimensions"); + } + + matlab_csc_t a_csc = {0}; + matlab_csc_t z_csc = {0}; + + if (extract_matlab_csc(mx_a, &a_csc) != 0) { + mexErrMsgIdAndTxt("BinSparse:InvalidMatrix", + "Failed to extract CSC data from A"); + } + + bool have_zeros = false; + if (mx_zeros) { + if (extract_matlab_csc(mx_zeros, &z_csc) != 0) { + mexErrMsgIdAndTxt("BinSparse:InvalidZeros", + "Failed to extract CSC data from Zeros"); + } + have_zeros = true; + } + + bsp_matrix_t csc_matrix; + if (have_zeros) { + build_csc_merged(&a_csc, &z_csc, &csc_matrix); + } else { + build_csc_from_a(&a_csc, &csc_matrix); + } + + bsp_matrix_t result = csc_matrix; + if (target_format != BSP_CSC) { + result = bsp_convert_matrix_allocator(csc_matrix, target_format, + bsp_matlab_allocator); + bsp_destroy_matrix_t(&csc_matrix); + + if (result.format != target_format) { + bsp_destroy_matrix_t(&result); + mexErrMsgIdAndTxt("BinSparse:ConversionError", + "Failed to convert matrix to requested format"); + } + } + + if (!matrix_uses_allocator(&result, bsp_matlab_allocator)) { + bsp_matrix_t copied; + bsp_error_t error = + bsp_matrix_copy_with_allocator(&result, &copied, bsp_matlab_allocator); + bsp_destroy_matrix_t(&result); + if (error != BSP_SUCCESS) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to allocate MATLAB-owned matrix"); + } + result = copied; + } + + plhs[0] = bsp_matrix_to_matlab_struct(&result); + bsp_destroy_matrix_t(&result); +} diff --git a/bindings/matlab/binsparse_from_ssmc.m b/bindings/matlab/binsparse_from_ssmc.m new file mode 100644 index 0000000..31bce1c --- /dev/null +++ b/bindings/matlab/binsparse_from_ssmc.m @@ -0,0 +1,32 @@ +function matrix = binsparse_from_ssmc (A, Zeros, format) +%BINSPARSE_FROM_SSMC convert SuiteSparse A+Zeros to a Binsparse matrix struct +% +% Usage: +% matrix = binsparse_from_ssmc (A) % default format +% matrix = binsparse_from_ssmc (A, format) +% matrix = binsparse_from_ssmc (A, Zeros) +% matrix = binsparse_from_ssmc (A, Zeros, format) +% +% A is a sparse or dense MATLAB matrix. Zeros is an optional sparse matrix +% of the same size as A holding the pattern of explicit zero entries (as in +% a SuiteSparse Matrix Collection Problem.Zeros); its entries are stored in +% the output with the value zero. format is an optional format string: +% 'COO' (the default), 'COOR', 'CSC', or 'CSR' for sparse input, or +% 'DVEC', 'DMAT', or 'DMATC' for dense input. The result is a Binsparse +% matrix struct suitable for binsparse_write. +% +% Example: +% A = sparse ([1 3], [1 2], [5 7], 3, 3) ; +% Zeros = sparse (2, 3, 1, 3, 3) ; +% matrix = binsparse_from_ssmc (A, Zeros, 'CSC') ; +% +% See also binsparse_write, binsparse_minimize_types, generate_bsp_from_ssmc. + +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +% This .m file provides the help text for the binsparse_from_ssmc MEX +% function. + +error('binsparse_from_ssmc mexFunction not found; compile with build_matlab_bindings first'); diff --git a/bindings/matlab/binsparse_minimize_types.c b/bindings/matlab/binsparse_minimize_types.c new file mode 100644 index 0000000..04ef621 --- /dev/null +++ b/bindings/matlab/binsparse_minimize_types.c @@ -0,0 +1,448 @@ +/* + * SPDX-FileCopyrightText: 2024 Binsparse Developers + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +/** + * binsparse_minimize_types.c - Minimize value/index types in a MATLAB + * Binsparse matrix struct. + * + * Usage in MATLAB/Octave: + * matrix = binsparse_minimize_types(matrix) + */ + +#include "mex.h" +#include +#include +#include + +#include "matlab_bsp_helpers.h" + +static bool array_is_signed_integer(bsp_type_t type) { + return type == BSP_INT8 || type == BSP_INT16 || type == BSP_INT32 || + type == BSP_INT64; +} + +static bool array_is_integer(bsp_type_t type) { + return type == BSP_UINT8 || type == BSP_UINT16 || type == BSP_UINT32 || + type == BSP_UINT64 || array_is_signed_integer(type); +} + +static void minimize_iso_values(bsp_matrix_t* matrix) { + if (matrix->is_iso || matrix->values.size == 0) { + return; + } + + bsp_array_t* values = &matrix->values; + bool all_equal = true; + + switch (values->type) { + case BSP_UINT8: + case BSP_UINT16: + case BSP_UINT32: + case BSP_UINT64: { + uint64_t first = 0; + bsp_array_read((*values), 0, first); + for (size_t i = 1; i < values->size; i++) { + uint64_t current = 0; + bsp_array_read((*values), i, current); + if (current != first) { + all_equal = false; + break; + } + } + if (!all_equal) { + return; + } + + bsp_type_t out_type = values->type; + if (first == 1) { + out_type = BSP_BINT8; + } + + bsp_array_t new_values; + bsp_error_t error = bsp_construct_array_t_allocator( + &new_values, 1, out_type, bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + return; + } + + if (out_type == BSP_BINT8) { + int8_t one = 1; + bsp_array_write(new_values, 0, one); + } else { + bsp_array_write(new_values, 0, first); + } + + bsp_destroy_array_t(values); + *values = new_values; + matrix->is_iso = true; + return; + } + case BSP_INT8: + case BSP_INT16: + case BSP_INT32: + case BSP_INT64: + case BSP_BINT8: { + int64_t first = 0; + bsp_array_read((*values), 0, first); + for (size_t i = 1; i < values->size; i++) { + int64_t current = 0; + bsp_array_read((*values), i, current); + if (current != first) { + all_equal = false; + break; + } + } + if (!all_equal) { + return; + } + + bsp_type_t out_type = values->type; + if (first == 1) { + out_type = BSP_BINT8; + } + + bsp_array_t new_values; + bsp_error_t error = bsp_construct_array_t_allocator( + &new_values, 1, out_type, bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + return; + } + + if (out_type == BSP_BINT8) { + int8_t one = 1; + bsp_array_write(new_values, 0, one); + } else { + bsp_array_write(new_values, 0, first); + } + + bsp_destroy_array_t(values); + *values = new_values; + matrix->is_iso = true; + return; + } + case BSP_FLOAT32: { + float first = 0.0f; + bsp_array_read((*values), 0, first); + for (size_t i = 1; i < values->size; i++) { + float current = 0.0f; + bsp_array_read((*values), i, current); + if (current != first) { + all_equal = false; + break; + } + } + if (!all_equal) { + return; + } + + bsp_array_t new_values; + bsp_error_t error = bsp_construct_array_t_allocator( + &new_values, 1, values->type, bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + return; + } + + bsp_array_write(new_values, 0, first); + bsp_destroy_array_t(values); + *values = new_values; + matrix->is_iso = true; + return; + } + case BSP_FLOAT64: { + double first = 0.0; + bsp_array_read((*values), 0, first); + for (size_t i = 1; i < values->size; i++) { + double current = 0.0; + bsp_array_read((*values), i, current); + if (current != first) { + all_equal = false; + break; + } + } + if (!all_equal) { + return; + } + + bsp_array_t new_values; + bsp_error_t error = bsp_construct_array_t_allocator( + &new_values, 1, values->type, bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + return; + } + + bsp_array_write(new_values, 0, first); + bsp_destroy_array_t(values); + *values = new_values; + matrix->is_iso = true; + return; + } + case BSP_COMPLEX_FLOAT32: { + float _Complex first = 0.0f + 0.0f * I; + bsp_array_read((*values), 0, first); + for (size_t i = 1; i < values->size; i++) { + float _Complex current = 0.0f + 0.0f * I; + bsp_array_read((*values), i, current); + if (current != first) { + all_equal = false; + break; + } + } + if (!all_equal) { + return; + } + + bsp_array_t new_values; + bsp_error_t error = bsp_construct_array_t_allocator( + &new_values, 1, values->type, bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + return; + } + + bsp_array_write(new_values, 0, first); + bsp_destroy_array_t(values); + *values = new_values; + matrix->is_iso = true; + return; + } + case BSP_COMPLEX_FLOAT64: { + double _Complex first = 0.0 + 0.0 * I; + bsp_array_read((*values), 0, first); + for (size_t i = 1; i < values->size; i++) { + double _Complex current = 0.0 + 0.0 * I; + bsp_array_read((*values), i, current); + if (current != first) { + all_equal = false; + break; + } + } + if (!all_equal) { + return; + } + + bsp_array_t new_values; + bsp_error_t error = bsp_construct_array_t_allocator( + &new_values, 1, values->type, bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + return; + } + + bsp_array_write(new_values, 0, first); + bsp_destroy_array_t(values); + *values = new_values; + matrix->is_iso = true; + return; + } + default: + return; + } +} + +static void minimize_int64_values(bsp_matrix_t* matrix) { + int64_t* values = (int64_t*) matrix->values.data; + + int64_t min_value = values[0]; + int64_t max_value = values[0]; + + for (size_t i = 1; i < matrix->values.size; i++) { + if (values[i] > max_value) { + max_value = values[i]; + } + if (values[i] < min_value) { + min_value = values[i]; + } + } + + bsp_type_t value_type = BSP_INT64; + if (min_value >= 0) { + value_type = bsp_pick_integer_type((size_t) max_value); + } else { + if (max_value <= (int64_t) INT8_MAX && min_value >= (int64_t) INT8_MIN) { + value_type = BSP_INT8; + } else if (max_value <= (int64_t) INT16_MAX && + min_value >= (int64_t) INT16_MIN) { + value_type = BSP_INT16; + } else if (max_value <= (int64_t) INT32_MAX && + min_value >= (int64_t) INT32_MIN) { + value_type = BSP_INT32; + } else { + value_type = BSP_INT64; + } + } + + if (value_type == matrix->values.type) { + return; + } + + bsp_array_t new_values; + bsp_error_t error = bsp_construct_array_t_allocator( + &new_values, matrix->values.size, value_type, bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + return; + } + + for (size_t i = 0; i < matrix->values.size; i++) { + int64_t value = values[i]; + bsp_array_write(new_values, i, value); + } + + bsp_destroy_array_t(&matrix->values); + matrix->values = new_values; +} + +static void minimize_values_matlab(bsp_matrix_t* matrix) { + if (matrix->values.size == 0) { + return; + } + + if (matrix->values.type == BSP_FLOAT64) { + double* values = (double*) matrix->values.data; + bool float32_representable = true; + for (size_t i = 0; i < matrix->values.size; i++) { + if (((float) values[i]) != values[i]) { + float32_representable = false; + break; + } + } + + if (float32_representable) { + bsp_array_t new_values; + bsp_error_t error = bsp_construct_array_t_allocator( + &new_values, matrix->values.size, BSP_FLOAT32, bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + return; + } + + float* n_values = (float*) new_values.data; + for (size_t i = 0; i < matrix->values.size; i++) { + n_values[i] = (float) values[i]; + } + + bsp_destroy_array_t(&matrix->values); + matrix->values = new_values; + } + } else if (matrix->values.type == BSP_INT64) { + minimize_int64_values(matrix); + } else if (matrix->values.type == BSP_COMPLEX_FLOAT64) { + bool float32_representable = true; + double _Complex* values = (double _Complex*) matrix->values.data; + + for (size_t i = 0; i < matrix->values.size; i++) { + if (((float _Complex) values[i]) != values[i]) { + float32_representable = false; + break; + } + } + + if (float32_representable) { + bsp_array_t new_values; + bsp_error_t error = bsp_construct_array_t_allocator( + &new_values, matrix->values.size, BSP_COMPLEX_FLOAT32, + bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + return; + } + + float _Complex* n_values = (float _Complex*) new_values.data; + for (size_t i = 0; i < matrix->values.size; i++) { + n_values[i] = (float _Complex) values[i]; + } + + bsp_destroy_array_t(&matrix->values); + matrix->values = new_values; + } + } +} + +static void minimize_index_array(bsp_array_t* array, const char* name) { + if (array->size == 0) { + return; + } + + if (!array_is_integer(array->type)) { + mexErrMsgIdAndTxt("BinSparse:InvalidIndexType", + "%s must be an integer array", name); + } + + size_t max_value = 0; + + if (array_is_signed_integer(array->type)) { + for (size_t i = 0; i < array->size; i++) { + int64_t value = 0; + bsp_array_read((*array), i, value); + if (value < 0) { + mexErrMsgIdAndTxt("BinSparse:InvalidIndexValue", + "%s contains negative values", name); + } + if ((uint64_t) value > max_value) { + max_value = (size_t) value; + } + } + } else { + for (size_t i = 0; i < array->size; i++) { + uint64_t value = 0; + bsp_array_read((*array), i, value); + if (value > max_value) { + max_value = (size_t) value; + } + } + } + + bsp_type_t target_type = bsp_pick_integer_type(max_value); + if (array->type == target_type) { + return; + } + + bsp_array_t new_array; + bsp_error_t error = bsp_construct_array_t_allocator( + &new_array, array->size, target_type, bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", "Failed to allocate %s array", + name); + } + + for (size_t i = 0; i < array->size; i++) { + uint64_t value = 0; + bsp_array_read((*array), i, value); + bsp_array_write(new_array, i, value); + } + + bsp_destroy_array_t(array); + *array = new_array; +} + +void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { + if (nrhs != 1) { + mexErrMsgIdAndTxt("BinSparse:InvalidArgs", + "Usage: matrix = binsparse_minimize_types(matrix)"); + } + + if (nlhs > 1) { + mexErrMsgIdAndTxt("BinSparse:TooManyOutputs", "Too many output arguments"); + } + + if (!mxIsStruct(prhs[0])) { + mexErrMsgIdAndTxt("BinSparse:InvalidMatrix", + "Input must be a Binsparse matrix struct"); + } + + bsp_matrix_t matrix; + bsp_error_t error = matlab_struct_to_bsp_matrix_allocator( + prhs[0], &matrix, bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + mexErrMsgIdAndTxt("BinSparse:ConversionError", + "Failed to convert MATLAB struct to matrix: %s", + bsp_get_error_string(error)); + } + + minimize_values_matlab(&matrix); + minimize_index_array(&matrix.indices_0, "indices_0"); + minimize_index_array(&matrix.indices_1, "indices_1"); + minimize_index_array(&matrix.pointers_to_1, "pointers_to_1"); + minimize_iso_values(&matrix); + + plhs[0] = bsp_matrix_to_matlab_struct(&matrix); + bsp_destroy_matrix_t(&matrix); +} diff --git a/bindings/matlab/binsparse_read.c b/bindings/matlab/binsparse_read.c new file mode 100644 index 0000000..6740797 --- /dev/null +++ b/bindings/matlab/binsparse_read.c @@ -0,0 +1,107 @@ +/* + * SPDX-FileCopyrightText: 2024 Binsparse Developers + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +/** + * binsparse_read.c - Read Binsparse matrices into MATLAB + * + * This MEX function reads Binsparse format matrices and returns them + * as MATLAB structs compatible with bsp_matrix_create. + * + * Usage in MATLAB/Octave: + * matrix = binsparse_read(filename) + * matrix = binsparse_read(filename, group) + */ + +#include "mex.h" +#include +#include +#include + +#include "matlab_bsp_helpers.h" + +// Keep this MEX function loaded for the whole MATLAB session: the HDF5 +// library used by libbinsparse installs process-wide state that is not safe +// to tear down when a MEX file is cleared. +static void lock_mex_module(void) { + if (!mexIsLocked()) { + mexLock(); + } +} + +/** + * Main MEX function entry point + */ +void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { + lock_mex_module(); + + char* filename = NULL; + char* group = NULL; + bsp_matrix_t matrix; + bsp_error_t error; + + // Check input arguments + if (nrhs < 1 || nrhs > 2) { + mexErrMsgIdAndTxt("BinSparse:InvalidArgs", + "Usage: matrix = binsparse_read(filename [, group])"); + } + + if (nlhs > 1) { + mexErrMsgIdAndTxt("BinSparse:TooManyOutputs", "Too many output arguments"); + } + + // Get filename + if (!mxIsChar(prhs[0])) { + mexErrMsgIdAndTxt("BinSparse:InvalidFilename", "Filename must be a string"); + } + + filename = mxArrayToString(prhs[0]); + if (!filename) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to convert filename string"); + } + + // Get optional group name + if (nrhs == 2) { + if (!mxIsChar(prhs[1])) { + mxFree(filename); + mexErrMsgIdAndTxt("BinSparse:InvalidGroup", + "Group name must be a string"); + } + + group = mxArrayToString(prhs[1]); + if (!group) { + mxFree(filename); + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to convert group string"); + } + } + + // Read the matrix using Binsparse + error = + bsp_read_matrix_allocator(&matrix, filename, group, bsp_matlab_allocator); + + if (error != BSP_SUCCESS) { + // Clean up + if (filename) + mxFree(filename); + if (group) + mxFree(group); + + mexErrMsgIdAndTxt("BinSparse:ReadError", "Failed to read matrix: %s", + bsp_get_error_string(error)); + } + + // Convert to MATLAB struct + plhs[0] = bsp_matrix_to_matlab_struct(&matrix); + + // Clean up + bsp_destroy_matrix_t(&matrix); + + if (filename) + mxFree(filename); + if (group) + mxFree(group); +} diff --git a/bindings/matlab/binsparse_read.m b/bindings/matlab/binsparse_read.m new file mode 100644 index 0000000..aba9650 --- /dev/null +++ b/bindings/matlab/binsparse_read.m @@ -0,0 +1,29 @@ +function matrix = binsparse_read (filename, group) +%BINSPARSE_READ read a matrix from a Binsparse HDF5 file +% +% Usage: +% matrix = binsparse_read (filename) % read from the file root +% matrix = binsparse_read (filename, group) % read from an HDF5 group +% +% filename is the path of a Binsparse HDF5 file (typically *.bsp.h5), and +% group optionally names an HDF5 group within the file that holds a +% Binsparse matrix. The result is a raw Binsparse matrix struct with the +% fields values, indices_0, indices_1, pointers_to_1, nrows, ncols, nnz, +% is_iso, format, and structure, as documented in bsp_matrix_create. +% Value and index arrays keep the types stored in the file; indices are +% 0-based. Use convert_to_problem_struct to build a MATLAB-native +% SuiteSparse Problem struct from Binsparse data. +% +% Example: +% matrix = binsparse_read ('west0067.bsp.h5') ; +% b = binsparse_read ('west0067.bsp.h5', 'b') ; +% +% See also binsparse_write, bsp_matrix_create, convert_to_problem_struct. + +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +% This .m file provides the help text for the binsparse_read MEX function. + +error ('binsparse_read mexFunction not found; compile with build_matlab_bindings first') ; diff --git a/bindings/matlab/binsparse_write.c b/bindings/matlab/binsparse_write.c new file mode 100644 index 0000000..e3d839d --- /dev/null +++ b/bindings/matlab/binsparse_write.c @@ -0,0 +1,160 @@ +/* + * SPDX-FileCopyrightText: 2024 Binsparse Developers + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +/** + * binsparse_write.c - Write MATLAB structs to Binsparse format + * + * This MEX function writes MATLAB structs (compatible with bsp_matrix_create) + * to Binsparse format files. + * + * Usage in MATLAB/Octave: + * binsparse_write(filename, matrix) + * binsparse_write(filename, matrix, group) + * binsparse_write(filename, matrix, group, json_string) + * binsparse_write(filename, matrix, group, json_string, compression_level) + */ + +#include "mex.h" +#include +#include +#include + +#include "matlab_bsp_helpers.h" + +// Keep this MEX function loaded for the whole MATLAB session: the HDF5 +// library used by libbinsparse installs process-wide state that is not safe +// to tear down when a MEX file is cleared. +static void lock_mex_module(void) { + if (!mexIsLocked()) { + mexLock(); + } +} + +/** + * Main MEX function entry point + */ +void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { + lock_mex_module(); + + char* filename = NULL; + char* group = NULL; + char* json_string = NULL; + int compression_level = 1; // Default compression + bsp_matrix_t matrix; + bsp_error_t error; + + // Check input arguments + if (nrhs < 2 || nrhs > 5) { + mexErrMsgIdAndTxt("BinSparse:InvalidArgs", + "Usage: binsparse_write(filename, matrix [, group [, " + "json_string [, compression_level]]])"); + } + + if (nlhs > 0) { + mexErrMsgIdAndTxt("BinSparse:TooManyOutputs", + "No output arguments expected"); + } + + // Get filename + if (!mxIsChar(prhs[0])) { + mexErrMsgIdAndTxt("BinSparse:InvalidFilename", "Filename must be a string"); + } + + filename = mxArrayToString(prhs[0]); + if (!filename) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to convert filename string"); + } + + // Get matrix struct + if (!mxIsStruct(prhs[1])) { + mxFree(filename); + mexErrMsgIdAndTxt("BinSparse:InvalidMatrix", "Matrix must be a struct"); + } + + // Convert MATLAB struct to bsp_matrix_t + error = matlab_struct_to_bsp_matrix_allocator(prhs[1], &matrix, + bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + mxFree(filename); + mexErrMsgIdAndTxt("BinSparse:ConversionError", + "Failed to convert MATLAB struct to matrix: %s", + bsp_get_error_string(error)); + } + + // Get optional group name + if (nrhs >= 3 && !mxIsEmpty(prhs[2])) { + if (!mxIsChar(prhs[2])) { + bsp_destroy_matrix_t(&matrix); + mxFree(filename); + mexErrMsgIdAndTxt("BinSparse:InvalidGroup", + "Group name must be a string"); + } + + group = mxArrayToString(prhs[2]); + if (!group) { + bsp_destroy_matrix_t(&matrix); + mxFree(filename); + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to convert group string"); + } + } + + // Get optional JSON string + if (nrhs >= 4 && !mxIsEmpty(prhs[3])) { + if (!mxIsChar(prhs[3])) { + bsp_destroy_matrix_t(&matrix); + if (group) + mxFree(group); + mxFree(filename); + mexErrMsgIdAndTxt("BinSparse:InvalidJSON", "JSON must be a string"); + } + + json_string = mxArrayToString(prhs[3]); + if (!json_string) { + bsp_destroy_matrix_t(&matrix); + if (group) + mxFree(group); + mxFree(filename); + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to convert JSON string"); + } + } + + // Get optional compression level + if (nrhs >= 5 && !mxIsEmpty(prhs[4])) { + if (!mxIsNumeric(prhs[4]) || mxIsComplex(prhs[4]) || + mxGetNumberOfElements(prhs[4]) != 1) { + bsp_destroy_matrix_t(&matrix); + if (json_string) + mxFree(json_string); + if (group) + mxFree(group); + mxFree(filename); + mexErrMsgIdAndTxt("BinSparse:InvalidCompression", + "Compression level must be a scalar integer"); + } + + compression_level = (int) mxGetScalar(prhs[4]); + } + + // Write the matrix using Binsparse + error = + bsp_write_matrix(filename, matrix, group, json_string, compression_level); + + // Clean up + bsp_destroy_matrix_t(&matrix); + if (json_string) + mxFree(json_string); + if (group) + mxFree(group); + mxFree(filename); + + if (error != BSP_SUCCESS) { + mexErrMsgIdAndTxt("BinSparse:WriteError", "Failed to write matrix: %s", + bsp_get_error_string(error)); + } +} diff --git a/bindings/matlab/binsparse_write.m b/bindings/matlab/binsparse_write.m new file mode 100644 index 0000000..c494ff1 --- /dev/null +++ b/bindings/matlab/binsparse_write.m @@ -0,0 +1,32 @@ +function binsparse_write(filename, matrix, group, json_string, compression_level) +%BINSPARSE_WRITE write a matrix to a file in Binsparse HDF5 format +% +% Usage: +% binsparse_write (filename, matrix) +% binsparse_write (filename, matrix, group) +% binsparse_write (filename, matrix, group, json_string) +% binsparse_write (filename, matrix, group, json_string, compression_level) +% +% filename is the path of the Binsparse HDF5 file to create (typically +% *.bsp.h5), and matrix is a Binsparse matrix struct as returned by +% binsparse_read, binsparse_from_ssmc, or bsp_matrix_create. The optional +% group names an HDF5 group to write into ('' or [ ] writes to the file +% root). The optional json_string is a JSON object whose keys are merged +% into the Binsparse descriptor as user metadata. The optional +% compression_level selects gzip compression from 0 (none) to 9; the +% default is 1. +% +% Example: +% matrix = binsparse_from_ssmc (sparse ([1 2], [1 2], [3 4])) ; +% binsparse_write ('example.bsp.h5', matrix) ; +% binsparse_write ('example.bsp.h5', matrix, 'b', '{"role": "b"}', 9) ; +% +% See also binsparse_read, binsparse_from_ssmc, bsp_matrix_create, +% generate_bsp_from_ssmc. + +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% SPDX-License-Identifier: BSD-3-Clause + +% This .m file provides the help text for the binsparse_write MEX function. + +error ('binsparse_write mexFunction not found; compile with build_matlab_bindings first') ; diff --git a/bindings/matlab/binsparse_write_string_dataset.c b/bindings/matlab/binsparse_write_string_dataset.c new file mode 100644 index 0000000..098c2db --- /dev/null +++ b/bindings/matlab/binsparse_write_string_dataset.c @@ -0,0 +1,200 @@ +/* + * SPDX-FileCopyrightText: 2024 Binsparse Developers + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +/** + * binsparse_write_string_dataset.c - Write an HDF5 UTF-8 string dataset. + * + * Usage in MATLAB/Octave: + * binsparse_write_string_dataset(filename, dataset_name, value) + * + * The value may be a character vector or a cell array of character vectors. + * Strings are stored as variable-length UTF-8 HDF5 strings, matching the usual + * h5py representation for Python str values. + */ + +#include "mex.h" +#include +#include +#include +#include + +// Keep this MEX function loaded for the whole MATLAB session, and stop HDF5 +// from registering atexit handlers: both guard against crashes when MATLAB +// tears down MEX files that share HDF5 process-wide state. +static void lock_mex_module(void) { + if (!mexIsLocked()) { + H5dont_atexit(); + mexLock(); + } +} + +static char* get_required_string(const mxArray* value, const char* name) { + if (!mxIsChar(value)) { + mexErrMsgIdAndTxt("BinSparse:InvalidString", + "%s must be a character vector", name); + } + + char* string = mxArrayToString(value); + if (!string) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", "Failed to read %s", name); + } + return string; +} + +static char** collect_strings(const mxArray* value, size_t* count, + bool* scalar_dataset) { + char** strings = NULL; + + if (mxIsCell(value)) { + *count = mxGetNumberOfElements(value); + *scalar_dataset = false; + strings = (char**) mxCalloc(*count, sizeof(char*)); + for (size_t i = 0; i < *count; i++) { + const mxArray* cell = mxGetCell(value, i); + if (!cell || !mxIsChar(cell)) { + mexErrMsgIdAndTxt("BinSparse:InvalidValue", + "Cell values must be character vectors"); + } + strings[i] = mxArrayToString(cell); + if (!strings[i]) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to read string cell"); + } + } + } else if (mxIsChar(value)) { + *count = 1; + *scalar_dataset = true; + strings = (char**) mxCalloc(1, sizeof(char*)); + strings[0] = mxArrayToString(value); + if (!strings[0]) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", "Failed to read string"); + } + } else { + mexErrMsgIdAndTxt("BinSparse:InvalidValue", + "Value must be a character vector or cellstr"); + } + + return strings; +} + +static void free_strings(char** strings, size_t count) { + if (!strings) { + return; + } + for (size_t i = 0; i < count; i++) { + if (strings[i]) { + mxFree(strings[i]); + } + } + mxFree(strings); +} + +void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { + (void) plhs; + lock_mex_module(); + + if (nrhs != 3) { + mexErrMsgIdAndTxt("BinSparse:InvalidArgs", + "Usage: binsparse_write_string_dataset(filename, " + "dataset_name, value)"); + } + + if (nlhs > 0) { + mexErrMsgIdAndTxt("BinSparse:TooManyOutputs", + "No output arguments expected"); + } + + char* filename = get_required_string(prhs[0], "filename"); + char* dataset_name = get_required_string(prhs[1], "dataset_name"); + + size_t count = 0; + bool scalar_dataset = false; + char** strings = collect_strings(prhs[2], &count, &scalar_dataset); + + hid_t file = H5I_INVALID_HID; + hid_t type = H5I_INVALID_HID; + hid_t space = H5I_INVALID_HID; + hid_t dset = H5I_INVALID_HID; + const char* error_id = NULL; + const char* error_message = NULL; + + if (access(filename, F_OK) == 0) { + file = H5Fopen(filename, H5F_ACC_RDWR, H5P_DEFAULT); + } else { + file = H5Fcreate(filename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + } + if (file == H5I_INVALID_HID) { + error_id = "BinSparse:FileError"; + error_message = "Failed to open HDF5 file"; + goto cleanup; + } + + htri_t exists = H5Lexists(file, dataset_name, H5P_DEFAULT); + if (exists > 0) { + error_id = "BinSparse:DatasetExists"; + error_message = "String dataset name already exists"; + goto cleanup; + } else if (exists < 0) { + error_id = "BinSparse:HDF5Error"; + error_message = "Failed to check string dataset name"; + goto cleanup; + } + + type = H5Tcopy(H5T_C_S1); + if (type == H5I_INVALID_HID || H5Tset_size(type, H5T_VARIABLE) < 0 || + H5Tset_cset(type, H5T_CSET_UTF8) < 0) { + error_id = "BinSparse:HDF5Error"; + error_message = "Failed to create UTF-8 string datatype"; + goto cleanup; + } + + if (scalar_dataset) { + space = H5Screate(H5S_SCALAR); + } else { + hsize_t dims[1] = {(hsize_t) count}; + space = H5Screate_simple(1, dims, NULL); + } + if (space == H5I_INVALID_HID) { + error_id = "BinSparse:HDF5Error"; + error_message = "Failed to create string dataspace"; + goto cleanup; + } + + dset = H5Dcreate2(file, dataset_name, type, space, H5P_DEFAULT, H5P_DEFAULT, + H5P_DEFAULT); + if (dset == H5I_INVALID_HID) { + error_id = "BinSparse:HDF5Error"; + error_message = "Failed to create string dataset"; + goto cleanup; + } + + if (H5Dwrite(dset, type, H5S_ALL, H5S_ALL, H5P_DEFAULT, strings) < 0) { + error_id = "BinSparse:HDF5Error"; + error_message = "Failed to write string dataset"; + goto cleanup; + } + +cleanup: + if (dset != H5I_INVALID_HID) { + H5Dclose(dset); + } + if (space != H5I_INVALID_HID) { + H5Sclose(space); + } + if (type != H5I_INVALID_HID) { + H5Tclose(type); + } + if (file != H5I_INVALID_HID) { + H5Fclose(file); + } + free_strings(strings, count); + mxFree(dataset_name); + mxFree(filename); + + if (error_id) { + mexErrMsgIdAndTxt(error_id, "%s", error_message); + } +} diff --git a/bindings/matlab/bsp_matrix_create.m b/bindings/matlab/bsp_matrix_create.m new file mode 100644 index 0000000..dc6f609 --- /dev/null +++ b/bindings/matlab/bsp_matrix_create.m @@ -0,0 +1,71 @@ +function matrix = bsp_matrix_create(varargin) +% BSP_MATRIX_CREATE - Create a Binsparse matrix struct +% +% Creates a MATLAB struct analogous to the C bsp_matrix_t structure. +% The struct contains four native MATLAB arrays and metadata fields. +% +% Usage: +% matrix = bsp_matrix_create() % Empty matrix +% matrix = bsp_matrix_create(values, indices_0, indices_1, pointers_to_1, ... +% nrows, ncols, nnz, is_iso, format, structure) +% +% Fields: +% values - MATLAB array of matrix values +% indices_0 - MATLAB array of first dimension indices +% indices_1 - MATLAB array of second dimension indices +% pointers_to_1 - MATLAB array of pointers for compressed formats +% nrows - Number of rows (integer) +% ncols - Number of columns (integer) +% nnz - Number of non-zeros (integer) +% is_iso - Logical indicating if matrix has single value (logical) +% format - Matrix format string ('CSR', 'CSC', 'COO', etc.) +% structure - Matrix structure string ('general', 'symmetric', etc.) +% +% Example: +% % Create empty matrix +% matrix = bsp_matrix_create(); +% +% % Create COO matrix +% values = [1.0, 2.0, 3.0]; +% rows = [1, 2, 3]; +% cols = [1, 2, 3]; +% matrix = bsp_matrix_create(values, rows, cols, [], 3, 3, 3, false, 'COO', 'general'); + +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +if nargin == 0 + % Create empty/default matrix + matrix = struct(... + 'values', double([]), ... + 'indices_0', uint64([]), ... + 'indices_1', uint64([]), ... + 'pointers_to_1', uint64([]), ... + 'nrows', uint64(0), ... + 'ncols', uint64(0), ... + 'nnz', uint64(0), ... + 'is_iso', false, ... + 'format', '', ... + 'structure', 'general'); + +elseif nargin == 10 + % Create matrix with all fields specified + matrix = struct(... + 'values', varargin{1}, ... + 'indices_0', varargin{2}, ... + 'indices_1', varargin{3}, ... + 'pointers_to_1', varargin{4}, ... + 'nrows', uint64(varargin{5}), ... + 'ncols', uint64(varargin{6}), ... + 'nnz', uint64(varargin{7}), ... + 'is_iso', logical(varargin{8}), ... + 'format', char(varargin{9}), ... + 'structure', char(varargin{10})); + +else + error('bsp_matrix_create:InvalidArgs', ... + 'Expected 0 or 10 arguments, got %d', nargin); +end + +end diff --git a/bindings/matlab/bsp_matrix_info.m b/bindings/matlab/bsp_matrix_info.m new file mode 100644 index 0000000..c2f0ed6 --- /dev/null +++ b/bindings/matlab/bsp_matrix_info.m @@ -0,0 +1,51 @@ +function bsp_matrix_info(matrix) +% BSP_MATRIX_INFO - Display information about a Binsparse matrix struct +% +% Usage: +% bsp_matrix_info(matrix) +% +% Displays: +% - Matrix dimensions and number of non-zeros +% - Format and structure information +% - Array sizes and types for each field + +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +if ~isstruct(matrix) + error('bsp_matrix_info:InvalidInput', 'Input must be a struct'); +end + +% Display basic matrix information +fprintf('%lu x %lu matrix with %lu nnz\n', ... + matrix.nrows, matrix.ncols, matrix.nnz); +fprintf('%s format with %s structure\n', ... + matrix.format, matrix.structure); + +if matrix.is_iso + fprintf('ISO matrix (single value)\n'); +end + +% Display array information +if ~isempty(matrix.values) + fprintf('%lu values of type %s\n', ... + length(matrix.values), class(matrix.values)); +end + +if ~isempty(matrix.indices_0) + fprintf('%lu indices_0 of type %s\n', ... + length(matrix.indices_0), class(matrix.indices_0)); +end + +if ~isempty(matrix.indices_1) + fprintf('%lu indices_1 of type %s\n', ... + length(matrix.indices_1), class(matrix.indices_1)); +end + +if ~isempty(matrix.pointers_to_1) + fprintf('%lu pointers_to_1 of type %s\n', ... + length(matrix.pointers_to_1), class(matrix.pointers_to_1)); +end + +end diff --git a/bindings/matlab/build_matlab_bindings.m b/bindings/matlab/build_matlab_bindings.m new file mode 100644 index 0000000..34408ca --- /dev/null +++ b/bindings/matlab/build_matlab_bindings.m @@ -0,0 +1,190 @@ +function build_matlab_bindings(varargin) +% BUILD_MATLAB_BINDINGS - Build Binsparse MATLAB MEX functions +% +% This script provides a simple interface to build MATLAB bindings +% for the Binsparse library. It automatically detects include paths +% and sets up the compilation environment. +% +% Usage: +% build_matlab_bindings() % Build all available MEX functions +% build_matlab_bindings('verbose') % Build with verbose output +% build_matlab_bindings('clean') % Clean compiled MEX files +% +% Prerequisites: +% - MATLAB with working MEX compiler (run 'mex -setup' if needed) +% - Binsparse C library headers (in ../../include/) +% - Compiled Binsparse library (in ../../build/) + +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +% Parse input arguments +verbose = any(strcmpi(varargin, 'verbose')); +clean_only = any(strcmpi(varargin, 'clean')); + +fprintf('=== Binsparse MATLAB Bindings Build Script ===\n\n'); + +if clean_only + fprintf('Cleaning compiled MEX files...\n'); + clean_mex_files(); + return; +end + +% Check MEX compiler +if ~check_mex_compiler() + error('MEX compiler not properly configured. Run "mex -setup" first.'); +end + +% Find and validate paths +paths = get_build_paths(); +if verbose + fprintf('Build paths:\n'); + fprintf(' MATLAB dir: %s\n', paths.matlab_dir); + fprintf(' Include dir: %s\n', paths.include_dir); + if ~isempty(paths.hdf5_include_dir) + fprintf(' HDF5 include dir: %s\n', paths.hdf5_include_dir); + end + fprintf(' Root dir: %s\n', paths.binsparse_root); + fprintf('\n'); +end + +% Compile MEX functions +compile_mex_functions(paths, verbose); + +fprintf('\n=== Build Complete ===\n'); +fprintf('Run the test functions to verify the installation:\n'); +fprintf(' test_binsparse_read()\n'); +fprintf(' test_binsparse_write()\n'); +fprintf(' test_write_binsparse_from_matlab()\n\n'); + +end + +function success = check_mex_compiler() + % Check if MEX compiler is configured + try + % Try to get MEX configuration + cc = mex.getCompilerConfigurations('C'); + success = ~isempty(cc); + if success + fprintf('MEX compiler found: %s\n', cc(1).Name); + end + catch + success = false; + end +end + +function paths = get_build_paths() + % Get and validate build paths + paths.matlab_dir = pwd; + paths.binsparse_root = fullfile(paths.matlab_dir, '..', '..'); + paths.include_dir = fullfile(paths.binsparse_root, 'include'); + paths.hdf5_include_dir = hdf5_include_dir(); + + if ~exist(paths.include_dir, 'dir') + error('Binsparse include directory not found: %s\nEnsure you are running this script from the bindings/matlab directory.', paths.include_dir); + end + + % Check for main header file + main_header = fullfile(paths.include_dir, 'binsparse', 'binsparse.h'); + if ~exist(main_header, 'file') + error('Main Binsparse header not found: %s', main_header); + end +end + +function include_dir = hdf5_include_dir() + include_dir = ''; + if exist('/usr/include/hdf5/serial', 'dir') + include_dir = '/usr/include/hdf5/serial'; + elseif exist('/usr/include/hdf5', 'dir') + include_dir = '/usr/include/hdf5'; + end +end + +function compile_mex_functions(paths, verbose) + % Compile all MEX functions + + % List of MEX functions to compile + mex_files = {'binsparse_read.c', 'binsparse_write.c', ... + 'binsparse_from_ssmc.c', 'binsparse_minimize_types.c', ... + 'write_binsparse_from_matlab.c', ... + 'binsparse_write_string_dataset.c'}; + + fprintf('Compiling MEX functions...\n'); + failed_files = {}; + + for i = 1:length(mex_files) + mex_file = mex_files{i}; + if ~exist(mex_file, 'file') + warning('MEX source file not found: %s', mex_file); + continue; + end + + fprintf(' Compiling %s... ', mex_file); + + % Prepare MEX command with library linking. The MEX functions link + % against the shared Binsparse library, so embed an rpath to the + % build directory on platforms that support it. + lib_dir = fullfile(paths.binsparse_root, 'build'); + lib_path = fullfile(lib_dir, 'libbinsparse_dynamic.so'); + cjson_lib = fullfile(lib_dir, '_deps', 'cjson-build', 'libcjson.a'); + if (ismac) + rpath = '-rpath ' ; + elseif (isunix) + rpath = '-rpath=' ; + else + rpath = '' ; + end + if ~isempty(rpath) + rpath = sprintf (' -Wl,%s''''%s'''' ', rpath, lib_dir) ; + rpath = [' LDFLAGS=''$LDFLAGS -fPIC ' rpath ' '' '] ; + end + + include_flags = sprintf('-I%s', paths.include_dir); + if ~isempty(paths.hdf5_include_dir) + include_flags = sprintf('%s -I%s', include_flags, ... + paths.hdf5_include_dir); + end + mex_command = sprintf ('mex %s %s %s %s %s -lhdf5_serial', ... + include_flags, rpath, lib_path, mex_file, cjson_lib) ; + if verbose + mex_command = [mex_command, ' -v'] ; + fprintf('\n %s\n', mex_command); + end + + try + eval (mex_command) ; + fprintf('SUCCESS\n'); + catch me + fprintf('FAILED\n'); + fprintf(' Error: %s\n', me.message); + failed_files{end+1} = mex_file; + end + end + + if ~isempty(failed_files) + error('Failed to compile MEX file(s): %s', ... + strjoin(failed_files, ', ')); + end +end + +function clean_mex_files() + % Clean compiled MEX files + + % Get MEX file extension for current platform + mex_ext = mexext(); + + % Find and delete MEX files + mex_files = dir(['*.' mex_ext]); + + if isempty(mex_files) + fprintf('No MEX files found to clean.\n'); + return; + end + + fprintf('Removing %d MEX file(s):\n', length(mex_files)); + for i = 1:length(mex_files) + fprintf(' %s\n', mex_files(i).name); + delete(mex_files(i).name); + end +end diff --git a/bindings/matlab/build_octave_bindings.m b/bindings/matlab/build_octave_bindings.m new file mode 100644 index 0000000..5031ac1 --- /dev/null +++ b/bindings/matlab/build_octave_bindings.m @@ -0,0 +1,213 @@ +function build_octave_bindings(varargin) +% BUILD_OCTAVE_BINDINGS - Build Binsparse Octave MEX functions +% +% This script provides a simple interface to build Octave bindings +% for the Binsparse library using mkoctfile. It automatically detects +% include paths and sets up the compilation environment. +% +% Usage: +% build_octave_bindings() % Build all available MEX functions +% build_octave_bindings('verbose') % Build with verbose output +% build_octave_bindings('clean') % Clean compiled MEX files +% +% Prerequisites: +% - GNU Octave with mkoctfile (usually included with Octave) +% - Binsparse C library headers (in ../../include/) +% - C compiler (gcc recommended) +% +% Note: This script builds Octave-compatible MEX functions using mkoctfile +% instead of MATLAB's mex command. + +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +% Parse input arguments +verbose = any(strcmpi(varargin, 'verbose')); +clean_only = any(strcmpi(varargin, 'clean')); + +fprintf('=== Binsparse Octave Bindings Build Script ===\n\n'); + +if clean_only + fprintf('Cleaning compiled MEX files...\n'); + clean_mex_files(); + return; +end + +% Check if we're running in Octave +if ~is_octave() + warning('This script is designed for Octave. For MATLAB, use build_matlab_bindings.m'); +end + +% Check mkoctfile availability +if ~check_mkoctfile() + error('mkoctfile not found. Please ensure Octave is properly installed.'); +end + +% Find and validate paths +paths = get_build_paths(); +if verbose + fprintf('Build paths:\n'); + fprintf(' Current dir: %s\n', paths.current_dir); + fprintf(' Include dir: %s\n', paths.include_dir); + if ~isempty(paths.hdf5_include_dir) + fprintf(' HDF5 include dir: %s\n', paths.hdf5_include_dir); + end + fprintf(' Root dir: %s\n', paths.binsparse_root); + fprintf('\n'); +end + +% Compile MEX functions +compile_octave_functions(paths, verbose); + +fprintf('\n=== Build Complete ===\n'); +fprintf('Run the test functions to verify the installation:\n'); +fprintf(' test_binsparse_read()\n'); +fprintf(' test_binsparse_write()\n'); +fprintf(' test_write_binsparse_from_matlab()\n\n'); + +end + +function result = is_octave() + % Check if running in Octave + result = exist('OCTAVE_VERSION', 'builtin') ~= 0; +end + +function success = check_mkoctfile() + % Check if mkoctfile is available + try + [status, ~] = system('mkoctfile --version'); + success = (status == 0); + if success && nargout == 0 + fprintf('mkoctfile found and working\n'); + end + catch + success = false; + end +end + +function paths = get_build_paths() + % Get and validate build paths + paths.current_dir = pwd; + paths.binsparse_root = fullfile(paths.current_dir, '..', '..'); + paths.include_dir = fullfile(paths.binsparse_root, 'include'); + paths.hdf5_include_dir = hdf5_include_dir(); + + if ~exist(paths.include_dir, 'dir') + error('Binsparse include directory not found: %s\nEnsure you are running this script from the bindings/matlab directory.', paths.include_dir); + end + + % Check for main header file + main_header = fullfile(paths.include_dir, 'binsparse', 'binsparse.h'); + if ~exist(main_header, 'file') + error('Main Binsparse header not found: %s', main_header); + end +end + +function include_dir = hdf5_include_dir() + include_dir = ''; + if exist('/usr/include/hdf5/serial', 'dir') + include_dir = '/usr/include/hdf5/serial'; + elseif exist('/usr/include/hdf5', 'dir') + include_dir = '/usr/include/hdf5'; + end +end + +function compile_octave_functions(paths, verbose) + % Compile all MEX functions using mkoctfile + + % List of MEX functions to compile + mex_files = {'binsparse_read.c', 'binsparse_write.c', ... + 'binsparse_from_ssmc.c', 'binsparse_minimize_types.c', ... + 'write_binsparse_from_matlab.c', ... + 'binsparse_write_string_dataset.c'}; + + fprintf('Compiling MEX functions with mkoctfile...\n'); + failed_files = {}; + + for i = 1:length(mex_files) + mex_file = mex_files{i}; + if ~exist(mex_file, 'file') + warning('MEX source file not found: %s', mex_file); + continue; + end + + fprintf(' Compiling %s... ', mex_file); + + % Prepare mkoctfile command with library linking + include_flag = sprintf('-I%s', paths.include_dir); + if ~isempty(paths.hdf5_include_dir) + include_flag = sprintf('%s -I%s', include_flag, ... + paths.hdf5_include_dir); + end + lib_dir = fullfile(paths.binsparse_root, 'build'); + lib_path = fullfile(lib_dir, 'libbinsparse.a'); + cjson_lib_dir = fullfile(lib_dir, '_deps', 'cjson-build'); + + if verbose + cmd = sprintf('mkoctfile --mex --verbose -fPIC %s %s -Wl,--whole-archive %s -Wl,--no-whole-archive -L%s -lcjson -lhdf5_serial', ... + include_flag, mex_file, lib_path, cjson_lib_dir); + else + cmd = sprintf('mkoctfile --mex -fPIC %s %s -Wl,--whole-archive %s -Wl,--no-whole-archive -L%s -lcjson -lhdf5_serial', ... + include_flag, mex_file, lib_path, cjson_lib_dir); + end + + if verbose + fprintf('\n Command: %s\n', cmd); + end + + % Execute mkoctfile + [status, output] = system(cmd); + + if status == 0 + fprintf('SUCCESS\n'); + if verbose && ~isempty(output) + fprintf(' Output: %s\n', output); + end + else + fprintf('FAILED\n'); + fprintf(' Error output:\n%s\n', output); + failed_files{end+1} = mex_file; + end + end + + if ~isempty(failed_files) + error('Failed to compile MEX file(s): %s', ... + strjoin(failed_files, ', ')); + end +end + +function clean_mex_files() + % Clean compiled MEX files (Octave uses different extensions) + + % Octave MEX extensions vary by platform + if ispc + extensions = {'mexw32', 'mexw64'}; + elseif ismac + extensions = {'mexmaci64'}; + else + extensions = {'mexa64', 'mex'}; + end + + found_files = {}; + + % Find files with any of the MEX extensions + for i = 1:length(extensions) + ext = extensions{i}; + files = dir(['*.' ext]); + for j = 1:length(files) + found_files{end+1} = files(j).name; + end + end + + if isempty(found_files) + fprintf('No MEX files found to clean.\n'); + return; + end + + fprintf('Removing %d MEX file(s):\n', length(found_files)); + for i = 1:length(found_files) + fprintf(' %s\n', found_files{i}); + delete(found_files{i}); + end +end diff --git a/bindings/matlab/compile_octave.sh b/bindings/matlab/compile_octave.sh new file mode 100755 index 0000000..4e429f2 --- /dev/null +++ b/bindings/matlab/compile_octave.sh @@ -0,0 +1,197 @@ +#!/bin/bash + +# SPDX-FileCopyrightText: 2024 Binsparse Developers +# +# SPDX-License-Identifier: BSD-3-Clause + +# compile_octave.sh - Build Binsparse Octave MEX functions from command line +# +# This script compiles the Binsparse MEX functions for Octave using mkoctfile +# from the command line, without needing to start Octave first. +# +# Usage: +# ./compile_octave.sh +# ./compile_octave.sh --verbose +# ./compile_octave.sh --clean + +set -e # Exit on any error + +# Color output for better readability +RED='\033[0;31m' +GREEN='\033[0;32m' +YELLOW='\033[1;33m' +BLUE='\033[0;34m' +NC='\033[0m' # No Color + +# Function to print colored output +print_info() { + echo -e "${BLUE}[INFO]${NC} $1" +} + +print_success() { + echo -e "${GREEN}[SUCCESS]${NC} $1" +} + +print_warning() { + echo -e "${YELLOW}[WARNING]${NC} $1" +} + +print_error() { + echo -e "${RED}[ERROR]${NC} $1" +} + +# Parse command line arguments +VERBOSE=false +CLEAN=false + +for arg in "$@"; do + case $arg in + --verbose|-v) + VERBOSE=true + shift + ;; + --clean|-c) + CLEAN=true + shift + ;; + --help|-h) + echo "Usage: $0 [--verbose] [--clean] [--help]" + echo " --verbose, -v: Enable verbose output" + echo " --clean, -c: Clean compiled MEX files" + echo " --help, -h: Show this help message" + exit 0 + ;; + *) + print_error "Unknown option: $arg" + echo "Use --help for usage information" + exit 1 + ;; + esac +done + +print_info "Binsparse Octave MEX Compilation Script" +echo "========================================" + +# Check if mkoctfile is available +if ! command -v mkoctfile &> /dev/null; then + print_error "mkoctfile not found. Please install GNU Octave." + exit 1 +fi + +mkoctfile_version=$(mkoctfile --version | head -n1) +print_info "Found mkoctfile: $mkoctfile_version" + +# Get script directory +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +print_info "Working directory: $SCRIPT_DIR" + +# Find Binsparse include directory +BINSPARSE_ROOT="$(realpath "$SCRIPT_DIR/../..")" +INCLUDE_DIR="$BINSPARSE_ROOT/include" + +if [ ! -d "$INCLUDE_DIR" ]; then + print_error "Binsparse include directory not found: $INCLUDE_DIR" + print_error "Make sure you're running this script from the bindings/matlab directory" + exit 1 +fi + +if [ ! -f "$INCLUDE_DIR/binsparse/binsparse.h" ]; then + print_error "Main Binsparse header not found: $INCLUDE_DIR/binsparse/binsparse.h" + exit 1 +fi + +print_info "Using Binsparse include directory: $INCLUDE_DIR" + +HDF5_INCLUDE_DIR="" +if [ -d /usr/include/hdf5/serial ]; then + HDF5_INCLUDE_DIR="/usr/include/hdf5/serial" +elif [ -d /usr/include/hdf5 ]; then + HDF5_INCLUDE_DIR="/usr/include/hdf5" +fi + +if [ -n "$HDF5_INCLUDE_DIR" ]; then + print_info "Using HDF5 include directory: $HDF5_INCLUDE_DIR" +fi + +# Change to script directory +cd "$SCRIPT_DIR" + +# Clean function +clean_mex_files() { + print_info "Cleaning compiled MEX files..." + + # MEX file extensions for different platforms + extensions=("mex" "mexa64" "mexw32" "mexw64" "mexmaci64") + + found_files=0 + for ext in "${extensions[@]}"; do + if ls *.${ext} 1> /dev/null 2>&1; then + for file in *.${ext}; do + print_info "Removing $file" + rm "$file" + ((found_files++)) + done + fi + done + + if [ $found_files -eq 0 ]; then + print_info "No MEX files found to clean" + else + print_success "Cleaned $found_files MEX file(s)" + fi +} + +# Handle clean option +if [ "$CLEAN" = true ]; then + clean_mex_files + exit 0 +fi + +# List of MEX files to compile +MEX_FILES=("binsparse_read.c" "binsparse_write.c" "binsparse_from_ssmc.c" "binsparse_minimize_types.c" "write_binsparse_from_matlab.c" "binsparse_write_string_dataset.c") + +print_info "Compiling MEX functions..." + +# Compile each MEX file +for mex_file in "${MEX_FILES[@]}"; do + if [ ! -f "$mex_file" ]; then + print_warning "MEX source file not found: $mex_file" + continue + fi + + print_info "Compiling $mex_file..." + + # Build mkoctfile command with library linking + LIB_DIR="$BINSPARSE_ROOT/build" + LIB_PATH="$LIB_DIR/libbinsparse.a" + CJSON_LIB_DIR="$LIB_DIR/_deps/cjson-build" + + INCLUDE_FLAGS="-I\"$INCLUDE_DIR\"" + if [ -n "$HDF5_INCLUDE_DIR" ]; then + INCLUDE_FLAGS="$INCLUDE_FLAGS -I\"$HDF5_INCLUDE_DIR\"" + fi + + CMD="mkoctfile --mex -fPIC $INCLUDE_FLAGS $mex_file -Wl,--whole-archive \"$LIB_PATH\" -Wl,--no-whole-archive -L\"$CJSON_LIB_DIR\" -lcjson -lhdf5_serial" + + if [ "$VERBOSE" = true ]; then + CMD="$CMD --verbose" + print_info "Command: $CMD" + fi + + # Execute compilation + if eval $CMD; then + print_success "Successfully compiled $mex_file" + else + print_error "Failed to compile $mex_file" + exit 1 + fi +done + +print_success "All MEX functions compiled successfully!" +echo "" +print_info "To test the functions, start Octave and run:" +echo " test_binsparse_read()" +echo " test_binsparse_write()" +echo "" +print_info "Or test from command line:" +echo " octave --eval \"test_binsparse_read()\"" diff --git a/bindings/matlab/convert_to_problem_struct.m b/bindings/matlab/convert_to_problem_struct.m new file mode 100644 index 0000000..4779cbf --- /dev/null +++ b/bindings/matlab/convert_to_problem_struct.m @@ -0,0 +1,492 @@ +function Problem = convert_to_problem_struct(bsp_problem) +%CONVERT_TO_PROBLEM_STRUCT convert a Binsparse problem to an SSMC Problem +% +% Problem = convert_to_problem_struct(bsp_problem) +% +% bsp_problem is an in-memory representation of one SuiteSparse Matrix +% Collection problem. Its A, b, x, and aux numeric entries are raw structs +% returned by binsparse_read. The metadata field contains the user metadata +% from the root Binsparse JSON descriptor. Text entries may be char, string, +% or cellstr values returned by h5read. + +% SPDX-FileCopyrightText: 2026 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +validate_problem(bsp_problem); + +Problem = metadata_to_problem(bsp_problem.metadata); +[Problem.A, Zeros] = convert_matrix(bsp_problem.A, true); +if ~issparse(Problem.A) || any(size(Problem.A) == 0) + error('BinSparse:InvalidProblem', ... + 'The primary SSMC matrix must be sparse and nonempty'); +end +if nnz(Zeros) > 0 + Problem.Zeros = Zeros; +end + +if isfield(bsp_problem, 'b') + Problem.b = convert_component(bsp_problem.b, Problem); +end +if isfield(bsp_problem, 'x') + Problem.x = convert_component(bsp_problem.x, Problem); +end +if isfield(bsp_problem, 'aux') + Problem.aux = convert_aux(bsp_problem.aux, Problem); +end + +end + +function validate_problem(bsp_problem) +if ~isstruct(bsp_problem) || ~isscalar(bsp_problem) + error('BinSparse:InvalidProblem', ... + 'Binsparse problem must be a scalar struct'); +end +if ~isfield(bsp_problem, 'metadata') || ... + ~isstruct(bsp_problem.metadata) || ~isscalar(bsp_problem.metadata) + error('BinSparse:InvalidProblem', ... + 'Binsparse problem must contain scalar metadata'); +end +if ~isfield(bsp_problem, 'A') || ~is_bsp_matrix(bsp_problem.A) + error('BinSparse:InvalidProblem', ... + 'Binsparse problem must contain a raw primary matrix A'); +end +end + +function Problem = metadata_to_problem(metadata) +Problem = struct(); +names = fieldnames(metadata); +for k = 1:numel(names) + name = names{k}; + if strcmp(name, 'role') + continue; + end + value = metadata.(name); + if strcmp(name, 'notes') + value = text_rows(value); + elseif any(strcmp(name, {'name', 'title', 'date', 'author', 'ed', ... + 'editor', 'kind', 'group', 'type', ... + 'structure'})) + value = text_scalar(value, name); + end + Problem.(name) = value; +end + +required = {'name', 'title', 'id', 'date', 'author', 'ed', 'kind'}; +for k = 1:numel(required) + if ~isfield(Problem, required{k}) + error('BinSparse:InvalidMetadata', ... + 'Missing required SSMC metadata field: %s', required{k}); + end +end +end + +function value = convert_component(value, Problem) +if is_bsp_matrix(value) + value = convert_matrix(value, false); +elseif is_text(value) + use_cellstr = isfield(Problem, 'id') && Problem.id > 2776; + value = normalize_component_text(value, use_cellstr); +else + error('BinSparse:InvalidComponent', ... + 'Unsupported Binsparse problem component'); +end +end + +function aux = convert_aux(raw_aux, Problem) +if ~isstruct(raw_aux) || ~isscalar(raw_aux) + error('BinSparse:InvalidAux', 'aux must be a scalar struct'); +end + +aux = struct(); +names = fieldnames(raw_aux); +for k = 1:numel(names) + name = names{k}; + value = convert_component(raw_aux.(name), Problem); + tokens = regexp(name, '^(.*)_([0-9]+)$', 'tokens', 'once'); + if isempty(tokens) || isempty(tokens{1}) + if isfield(aux, name) + error('BinSparse:DuplicateComponent', ... + 'Duplicate auxiliary component: %s', name); + end + aux.(name) = value; + continue; + end + + base = tokens{1}; + index = str2double(tokens{2}); + if ~isfinite(index) || index < 1 || fix(index) ~= index + error('BinSparse:InvalidComponent', ... + 'Invalid auxiliary cell component: %s', name); + end + if isfield(aux, base) && ~iscell(aux.(base)) + error('BinSparse:DuplicateComponent', ... + 'Auxiliary component conflicts with cell sequence: %s', base); + end + if ~isfield(aux, base) + aux.(base) = cell(0, 1); + end + if numel(aux.(base)) >= index && ~isempty(aux.(base){index}) + error('BinSparse:DuplicateComponent', ... + 'Duplicate auxiliary cell component: %s', name); + end + aux.(base){index, 1} = value; +end +end + +function [matrix, Zeros] = convert_matrix(bsp, split_zeros) +validate_matrix_header(bsp); +format = upper(char(bsp.format)); +m = checked_size(bsp.nrows, 'nrows'); +n = checked_size(bsp.ncols, 'ncols'); +count = checked_size(bsp.nnz, 'nnz'); + +switch format + case 'DVEC' + values = matrix_values(bsp, count); + if n ~= 1 || count ~= m + error('BinSparse:InvalidMatrix', ... + 'DVEC dimensions do not match the stored value count'); + end + matrix = reshape(values, [m, 1]); + Zeros = sparse(m, n); + return; + + case {'DMAT', 'DMATR'} + expected = checked_product(m, n); + if count ~= expected + error('BinSparse:InvalidMatrix', ... + 'DMATR dimensions do not match the stored value count'); + end + values = matrix_values(bsp, count); + matrix = reshape(values, [n, m]).'; + Zeros = sparse(m, n); + return; + + case 'DMATC' + expected = checked_product(m, n); + if count ~= expected + error('BinSparse:InvalidMatrix', ... + 'DMATC dimensions do not match the stored value count'); + end + values = matrix_values(bsp, count); + matrix = reshape(values, [m, n]); + Zeros = sparse(m, n); + return; +end + +[rows, cols] = sparse_indices(bsp, format, m, n, count); +values = matrix_values(bsp, count); +[rows, cols, values] = expand_structure(bsp, rows, cols, values, m, n); + +is_zero = (real(values) == 0) & (imag(values) == 0); +if split_zeros + matrix = sparse(rows(~is_zero), cols(~is_zero), ... + values(~is_zero), m, n); + Zeros = sparse(rows(is_zero), cols(is_zero), 1, m, n); +else + matrix = sparse(rows, cols, values, m, n); + Zeros = sparse(m, n); +end +end + +function [rows, cols] = sparse_indices(bsp, format, m, n, count) +switch format + case {'COO', 'COOR'} + rows = checked_indices(bsp.indices_0, count, m, 'indices_0') + 1; + cols = checked_indices(bsp.indices_1, count, n, 'indices_1') + 1; + require_ordered_pairs(rows, cols, 'COOR'); + + case 'COOC' + cols = checked_indices(bsp.indices_0, count, n, 'indices_0') + 1; + rows = checked_indices(bsp.indices_1, count, m, 'indices_1') + 1; + require_ordered_pairs(cols, rows, 'COOC'); + + case 'CSR' + pointers = checked_pointers(bsp.pointers_to_1, m + 1, count); + cols = checked_indices(bsp.indices_1, count, n, 'indices_1') + 1; + rows = expand_major_indices((0:m-1).', diff(pointers)) + 1; + require_segment_order(cols, pointers, 'CSR'); + + case 'CSC' + pointers = checked_pointers(bsp.pointers_to_1, n + 1, count); + rows = checked_indices(bsp.indices_1, count, m, 'indices_1') + 1; + cols = expand_major_indices((0:n-1).', diff(pointers)) + 1; + require_segment_order(rows, pointers, 'CSC'); + + case 'DCSR' + major = checked_indices(bsp.indices_0, [], m, 'indices_0'); + require_strictly_increasing(major, 'DCSR major indices'); + pointers = checked_pointers(bsp.pointers_to_1, numel(major) + 1, count); + require_nonempty_segments(pointers, 'DCSR'); + cols = checked_indices(bsp.indices_1, count, n, 'indices_1') + 1; + rows = expand_major_indices(major, diff(pointers)) + 1; + require_segment_order(cols, pointers, 'DCSR'); + + case 'DCSC' + major = checked_indices(bsp.indices_0, [], n, 'indices_0'); + require_strictly_increasing(major, 'DCSC major indices'); + pointers = checked_pointers(bsp.pointers_to_1, numel(major) + 1, count); + require_nonempty_segments(pointers, 'DCSC'); + rows = checked_indices(bsp.indices_1, count, m, 'indices_1') + 1; + cols = expand_major_indices(major, diff(pointers)) + 1; + require_segment_order(rows, pointers, 'DCSC'); + + case 'CVEC' + if n ~= 1 + error('BinSparse:InvalidMatrix', ... + 'CVEC must have one MATLAB column'); + end + rows = checked_indices(bsp.indices_0, count, m, 'indices_0') + 1; + require_strictly_increasing(rows, 'CVEC indices'); + cols = ones(count, 1); + + otherwise + error('BinSparse:UnsupportedFormat', ... + 'Unsupported Binsparse matrix format: %s', format); +end +end + +function values = matrix_values(bsp, count) +values = bsp.values(:); +if logical(bsp.is_iso) + if count == 0 + if numel(values) > 1 + error('BinSparse:InvalidMatrix', ... + 'An empty ISO matrix has too many values'); + end + values = values([]); + elseif numel(values) ~= 1 + error('BinSparse:InvalidMatrix', ... + 'An ISO matrix must contain exactly one value'); + else + values = repmat(values, count, 1); + end +elseif numel(values) ~= count + error('BinSparse:InvalidMatrix', ... + 'Value array length does not match nnz'); +end +values = exact_double(values); +end + +function [rows, cols, values] = expand_structure(bsp, rows, cols, values, m, n) +structure = 'general'; +if isfield(bsp, 'structure') && ~isempty(bsp.structure) + structure = lower(char(bsp.structure)); +end +if strcmp(structure, 'general') + return; +end +if m ~= n + error('BinSparse:InvalidStructure', ... + 'A structured Binsparse matrix must be square'); +end + +is_lower = ends_with(structure, '_lower'); +is_upper = ends_with(structure, '_upper'); +if ~is_lower && ~is_upper + error('BinSparse:UnsupportedStructure', ... + 'Unsupported Binsparse matrix structure: %s', structure); +end +if (is_lower && any(rows < cols)) || (is_upper && any(rows > cols)) + error('BinSparse:InvalidStructure', ... + 'Stored entries do not match the declared matrix triangle'); +end + +off_diagonal = (rows ~= cols); +mirror_rows = cols(off_diagonal); +mirror_cols = rows(off_diagonal); +mirror_values = values(off_diagonal); +if starts_with(structure, 'hermitian_') + mirror_values = conj(mirror_values); +elseif starts_with(structure, 'skew_symmetric_') + mirror_values = -mirror_values; +elseif ~starts_with(structure, 'symmetric_') + error('BinSparse:UnsupportedStructure', ... + 'Unsupported Binsparse matrix structure: %s', structure); +end + +rows = [rows; mirror_rows]; +cols = [cols; mirror_cols]; +values = [values; mirror_values]; +end + +function validate_matrix_header(bsp) +if ~is_bsp_matrix(bsp) || ~isscalar(bsp) + error('BinSparse:InvalidMatrix', ... + 'Expected a scalar raw Binsparse matrix struct'); +end +if ~isscalar(bsp.is_iso) || ... + (~islogical(bsp.is_iso) && ~any(bsp.is_iso == [0 1])) + error('BinSparse:InvalidMatrix', 'is_iso must be a logical scalar'); +end +end + +function values = exact_double(values) +converted = double(values); +if isinteger(values) && ~isequal(cast(converted, class(values)), values) + error('BinSparse:InexactValue', ... + 'Integer values cannot be represented exactly in MATLAB sparse form'); +end +values = converted; +end + +function ok = is_bsp_matrix(value) +required = {'values', 'indices_0', 'indices_1', 'pointers_to_1', ... + 'nrows', 'ncols', 'nnz', 'is_iso', 'format'}; +ok = isstruct(value) && isscalar(value) && ... + all(isfield(value, required)); +end + +function value = checked_size(value, name) +if ~isnumeric(value) || ~isreal(value) || ~isscalar(value) + error('BinSparse:InvalidMatrix', '%s must be a real numeric scalar', name); +end +value = double(value); +if ~isfinite(value) || value < 0 || fix(value) ~= value || value > flintmax + error('BinSparse:InvalidMatrix', 'Invalid %s value', name); +end +end + +function value = checked_product(a, b) +value = a * b; +if ~isfinite(value) || value > flintmax || fix(value) ~= value + error('BinSparse:InvalidMatrix', 'Dense matrix size is too large'); +end +end + +function indices = checked_indices(indices, expected, limit, name) +if ~isnumeric(indices) || ~isreal(indices) + error('BinSparse:InvalidMatrix', '%s must be real numeric data', name); +end +indices = double(indices(:)); +if ~isempty(expected) && numel(indices) ~= expected + error('BinSparse:InvalidMatrix', ... + '%s length does not match nnz', name); +end +if any(~isfinite(indices)) || any(indices < 0) || ... + any(fix(indices) ~= indices) || any(indices >= limit) + error('BinSparse:InvalidMatrix', '%s contains an invalid index', name); +end +end + +function pointers = checked_pointers(pointers, expected, count) +if ~isnumeric(pointers) || ~isreal(pointers) + error('BinSparse:InvalidMatrix', ... + 'pointers_to_1 must be real numeric data'); +end +pointers = double(pointers(:)); +if numel(pointers) ~= expected || isempty(pointers) || ... + pointers(1) ~= 0 || pointers(end) ~= count || ... + any(~isfinite(pointers)) || any(fix(pointers) ~= pointers) || ... + any(diff(pointers) < 0) + error('BinSparse:InvalidMatrix', 'Invalid pointers_to_1 array'); +end +end + +function expanded = expand_major_indices(major, counts) +if numel(major) ~= numel(counts) + error('BinSparse:InvalidMatrix', 'Major index and pointer size mismatch'); +end +expanded = repelem(double(major(:)), double(counts(:))); +expanded = expanded(:); +end + +function require_ordered_pairs(first, second, format) +if numel(first) < 2 + return; +end +bad = diff(first) < 0 | ... + (diff(first) == 0 & diff(second) <= 0); +if any(bad) + error('BinSparse:InvalidMatrix', ... + '%s indices are not sorted and unique', format); +end +end + +function require_segment_order(indices, pointers, format) +for k = 1:numel(pointers)-1 + first = pointers(k) + 1; + last = pointers(k + 1); + if last > first && any(diff(indices(first:last)) <= 0) + error('BinSparse:InvalidMatrix', ... + '%s minor indices are not sorted and unique', format); + end +end +end + +function require_nonempty_segments(pointers, format) +if any(diff(pointers) <= 0) + error('BinSparse:InvalidMatrix', ... + '%s must not store empty major dimensions', format); +end +end + +function require_strictly_increasing(values, label) +if numel(values) > 1 && any(diff(values) <= 0) + error('BinSparse:InvalidMatrix', '%s are not sorted and unique', label); +end +end + +function value = text_scalar(value, name) +if isstring(value) && isscalar(value) + value = char(value); +elseif iscell(value) && isscalar(value) + value = char(value{1}); +elseif ischar(value) && (isrow(value) || isempty(value)) + value = char(value); +else + error('BinSparse:InvalidMetadata', ... + 'Metadata field %s must be scalar text', name); +end +end + +function value = text_rows(value) +if ischar(value) + return; +elseif isstring(value) + value = cellstr(value(:)); +elseif iscellstr(value) + value = value(:); +else + error('BinSparse:InvalidMetadata', 'notes must contain text'); +end +value = char(value); +end + +function value = normalize_component_text(value, use_cellstr) +if ischar(value) + if size(value, 1) > 1 + rows = cellstr(value); + else + rows = {value}; + end +elseif isstring(value) + rows = cellstr(value(:)); +elseif iscellstr(value) + rows = value(:); +else + error('BinSparse:InvalidComponent', 'Text component is invalid'); +end +for k = 1:numel(rows) + rows{k} = reshape(char(rows{k}), 1, []); +end +if use_cellstr + value = rows; +else + value = char(rows); +end +end + +function ok = is_text(value) +ok = ischar(value) || iscellstr(value) || isstring(value); +end + +function ok = starts_with(value, prefix) +ok = strncmp(value, prefix, numel(prefix)); +end + +function ok = ends_with(value, suffix) +ok = numel(value) >= numel(suffix) && ... + strcmp(value(end-numel(suffix)+1:end), suffix); +end diff --git a/bindings/matlab/generate_bsp_from_ssmc.m b/bindings/matlab/generate_bsp_from_ssmc.m new file mode 100644 index 0000000..a7ac027 --- /dev/null +++ b/bindings/matlab/generate_bsp_from_ssmc.m @@ -0,0 +1,246 @@ +function generate_bsp_from_ssmc(problem, output_filename, format, compression_level) +% GENERATE_BSP_FROM_SSMC - Generate a Binsparse file from an SSMC problem +% +% Usage: +% generate_bsp_from_ssmc(problem, output_filename) +% generate_bsp_from_ssmc(problem, output_filename, format, compression_level) +% +% Defaults: +% format = 'COO' +% compression_level = 0 + +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +if nargin < 2 + error('generate_bsp_from_ssmc:InvalidArgs', ... + 'Usage: generate_bsp_from_ssmc(problem, output_filename [, format [, compression_level]])'); +end + +if nargin < 3 || isempty(format) + format = 'COO'; +end + +if nargin < 4 || isempty(compression_level) + compression_level = 0; +end + +if ~isstruct(problem) + error('generate_bsp_from_ssmc:InvalidProblem', 'Problem must be a struct'); +end + +if isfield(problem, 'Problem') + P = problem.Problem; +else + P = problem; +end + +if ~isfield(P, 'A') + error('generate_bsp_from_ssmc:MissingMatrix', 'Problem.A is required'); +end + +% Primary matrix +A = P.A; +if issparse(A) + if isfield(P, 'Zeros') && ~isempty(P.Zeros) + Zeros = P.Zeros; + mat = binsparse_from_ssmc(A, Zeros, format); + else + mat = binsparse_from_ssmc(A, format); + end +else + mat = binsparse_from_ssmc(A, dense_format_for(A)); +end + +mat = binsparse_minimize_types(mat); +primary_metadata = metadata_json(P, 'A'); +binsparse_write(output_filename, mat, '', primary_metadata, compression_level); + +% Handle aux struct +if isfield(P, 'aux') && isstruct(P.aux) + aux_names = fieldnames(P.aux); + for i = 1:numel(aux_names) + name = aux_names{i}; + value = P.aux.(name); + handle_aux_entry(name, value, output_filename, format, compression_level); + end +end + +% Handle x and b like aux +if isfield(P, 'x') && ~isempty(P.x) + handle_aux_entry('x', P.x, output_filename, format, compression_level); +end + +if isfield(P, 'b') && ~isempty(P.b) + handle_aux_entry('b', P.b, output_filename, format, compression_level); +end + +end + +function handle_aux_entry(name, value, output_filename, format, compression_level) + if isempty(value) + return; + end + + if is_text_value(value) + write_string_dataset(output_filename, name, value); + return; + end + + if iscell(value) + len = numel(value); + for k = 1:len + handle_aux_entry(component_name(name, k, len), value{k}, ... + output_filename, format, compression_level); + end + return; + end + + if islogical(value) + value = double(value); + end + + if issparse(value) + bsp = binsparse_from_ssmc(value, format); + elseif isnumeric(value) + bsp = binsparse_from_ssmc(value, dense_format_for(value)); + else + fprintf('Note: aux entry "%s" has unsupported type and is ignored.\n', name); + return; + end + + bsp = binsparse_minimize_types(bsp); + binsparse_write(output_filename, bsp, name, entry_metadata_json(name), ... + compression_level); +end + +function fmt = dense_format_for(value) + if isvector(value) + fmt = 'DVEC'; + else + fmt = 'DMATC'; + end +end + +function name = component_name(prefix, index, count) + if count < 10 + pattern = '%s_%d'; + elseif count < 100 + if index < 10 + pattern = '%s_0%d'; + else + pattern = '%s_%d'; + end + else + if index < 10 + pattern = '%s_00%d'; + elseif index < 100 + pattern = '%s_0%d'; + else + pattern = '%s_%d'; + end + end + name = sprintf(pattern, prefix, index); +end + +function ok = is_text_value(value) + ok = ischar(value) || isstring(value) || iscellstr(value); +end + +function write_string_dataset(output_filename, name, value) + if exist('binsparse_write_string_dataset', 'file') ~= 3 + error('generate_bsp_from_ssmc:MissingStringWriter', ... + 'BSP text output requires binsparse_write_string_dataset on the path'); + end + + if isstring(value) + if isscalar(value) + value = char(value); + else + value = cellstr(value(:)); + end + elseif ischar(value) + if size(value, 1) > 1 + value = char_rows(value); + else + value = char(value); + end + elseif iscellstr(value) + value = value(:); + else + error('generate_bsp_from_ssmc:InvalidStringValue', ... + 'Text aux value must be char, string, or cellstr'); + end + + binsparse_write_string_dataset(output_filename, name, value); +end + +function json = metadata_json(P, role) + metadata_fields = {'name', 'title', 'id', 'date', 'author', 'ed', ... + 'editor', 'kind', 'notes', 'group', 'num_rows', ... + 'num_cols', 'nnz', 'pattern_symmetry', ... + 'numerical_symmetry', 'type', 'structure'}; + metadata = struct(); + metadata.role = role; + + for i = 1:numel(metadata_fields) + field = metadata_fields{i}; + if isfield(P, field) + [ok, value] = metadata_value(P.(field)); + if ok + metadata.(field) = value; + end + end + end + + json = metadata_payload_json(metadata); +end + +function json = entry_metadata_json(role) + metadata = struct(); + metadata.role = role; + json = metadata_payload_json(metadata); +end + +function json = metadata_payload_json(metadata) + payload = struct(); + payload.metadata = metadata; + json = encode_json_or_empty(payload); +end + +function [ok, value] = metadata_value(value) + ok = false; + if ischar(value) + if size(value, 1) > 1 + value = char_rows(value); + else + value = char(value); + end + ok = true; + elseif isstring(value) + if isscalar(value) + value = char(value); + else + value = cellstr(value(:)); + end + ok = true; + elseif (islogical(value) || isnumeric(value)) && numel(value) <= 64 + ok = true; + end +end + +function rows = char_rows(value) + rows = cell(size(value, 1), 1); + for k = 1:size(value, 1) + rows{k} = value(k, :); + end +end + +function json = encode_json_or_empty(value) + try + json = jsonencode(value); + catch + json = []; + end +end diff --git a/bindings/matlab/matlab_bsp_helpers.h b/bindings/matlab/matlab_bsp_helpers.h new file mode 100644 index 0000000..02c607f --- /dev/null +++ b/bindings/matlab/matlab_bsp_helpers.h @@ -0,0 +1,484 @@ +/* SPDX-FileCopyrightText: 2024 Binsparse Developers + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#ifndef MATLAB_BSP_HELPERS_H +#define MATLAB_BSP_HELPERS_H + +#include "mex.h" +#include +#include +#include + +typedef struct { + double* values; + double* imag_values; + mwIndex* rowind; + mwIndex* colptr; + size_t nrows; + size_t ncols; + size_t nnz; + bool has_values; + bool is_complex; +} matlab_csc_t; + +static const bsp_allocator_t bsp_matlab_allocator = {.malloc = mxMalloc, + .free = mxFree}; + +static inline int extract_matlab_csc(const mxArray* mx_matrix, + matlab_csc_t* csc_matrix) { + if (!mx_matrix || !csc_matrix) { + return -1; + } + + if (!mxIsSparse(mx_matrix)) { + return -1; + } + + if (!mxIsDouble(mx_matrix) && !mxIsLogical(mx_matrix)) { + return -1; + } + + if (mxIsLogical(mx_matrix) && mxIsComplex(mx_matrix)) { + return -1; + } + + csc_matrix->nrows = mxGetM(mx_matrix); + csc_matrix->ncols = mxGetN(mx_matrix); + + csc_matrix->has_values = mxIsDouble(mx_matrix); + csc_matrix->is_complex = mxIsComplex(mx_matrix); + csc_matrix->values = csc_matrix->has_values ? mxGetPr(mx_matrix) : NULL; + csc_matrix->imag_values = csc_matrix->is_complex ? mxGetPi(mx_matrix) : NULL; + csc_matrix->rowind = mxGetIr(mx_matrix); + csc_matrix->colptr = mxGetJc(mx_matrix); + + if (!csc_matrix->rowind || !csc_matrix->colptr) { + return -1; + } + + if (csc_matrix->ncols > 0) { + csc_matrix->nnz = csc_matrix->colptr[csc_matrix->ncols]; + } else { + csc_matrix->nnz = 0; + } + + if (csc_matrix->nnz > 0 && csc_matrix->has_values && !csc_matrix->values) { + return -1; + } + + return 0; +} + +static inline bsp_error_t +matlab_to_bsp_array_allocator(const mxArray* mx_array, bsp_array_t* array, + bsp_allocator_t allocator) { + bool is_empty = mxIsEmpty(mx_array); + + size_t size = mxGetNumberOfElements(mx_array); + mxClassID class_id = mxGetClassID(mx_array); + bool is_complex = mxIsComplex(mx_array); + + bsp_type_t bsp_type; + size_t element_size; + + if (is_complex) { + if (class_id == mxDOUBLE_CLASS) { + bsp_type = BSP_COMPLEX_FLOAT64; + element_size = sizeof(double _Complex); + } else if (class_id == mxSINGLE_CLASS) { + bsp_type = BSP_COMPLEX_FLOAT32; + element_size = sizeof(float _Complex); + } else { + return BSP_INVALID_TYPE; + } + } else { + switch (class_id) { + case mxDOUBLE_CLASS: + bsp_type = BSP_FLOAT64; + element_size = sizeof(double); + break; + case mxSINGLE_CLASS: + bsp_type = BSP_FLOAT32; + element_size = sizeof(float); + break; + case mxUINT64_CLASS: + bsp_type = BSP_UINT64; + element_size = sizeof(uint64_t); + break; + case mxUINT32_CLASS: + bsp_type = BSP_UINT32; + element_size = sizeof(uint32_t); + break; + case mxUINT16_CLASS: + bsp_type = BSP_UINT16; + element_size = sizeof(uint16_t); + break; + case mxUINT8_CLASS: + bsp_type = BSP_UINT8; + element_size = sizeof(uint8_t); + break; + case mxINT64_CLASS: + bsp_type = BSP_INT64; + element_size = sizeof(int64_t); + break; + case mxINT32_CLASS: + bsp_type = BSP_INT32; + element_size = sizeof(int32_t); + break; + case mxINT16_CLASS: + bsp_type = BSP_INT16; + element_size = sizeof(int16_t); + break; + case mxINT8_CLASS: + bsp_type = BSP_INT8; + element_size = sizeof(int8_t); + break; + default: + return BSP_INVALID_TYPE; + } + } + + if (is_empty) { + array->data = NULL; + array->size = 0; + array->type = bsp_type; + array->allocator = allocator; + return BSP_SUCCESS; + } + + bsp_error_t error = + bsp_construct_array_t_allocator(array, size, bsp_type, allocator); + if (error != BSP_SUCCESS) { + return error; + } + + if (is_complex) { + if (class_id == mxDOUBLE_CLASS) { + double* real_data = mxGetPr(mx_array); + double* imag_data = mxGetPi(mx_array); + double _Complex* out_data = (double _Complex*) array->data; + for (size_t i = 0; i < size; i++) { + double imag = imag_data ? imag_data[i] : 0.0; + out_data[i] = real_data[i] + imag * I; + } + } else { + float* real_data = (float*) mxGetData(mx_array); + float* imag_data = (float*) mxGetImagData(mx_array); + float _Complex* out_data = (float _Complex*) array->data; + for (size_t i = 0; i < size; i++) { + float imag = imag_data ? imag_data[i] : 0.0f; + out_data[i] = real_data[i] + imag * I; + } + } + } else { + memcpy(array->data, mxGetData(mx_array), size * element_size); + } + + return BSP_SUCCESS; +} + +static inline bsp_error_t matlab_struct_to_bsp_matrix_allocator( + const mxArray* mx_struct, bsp_matrix_t* matrix, bsp_allocator_t allocator) { + bsp_construct_default_matrix_t_allocator(matrix, allocator); + + mxArray* values_field = mxGetField(mx_struct, 0, "values"); + mxArray* indices_0_field = mxGetField(mx_struct, 0, "indices_0"); + mxArray* indices_1_field = mxGetField(mx_struct, 0, "indices_1"); + mxArray* pointers_to_1_field = mxGetField(mx_struct, 0, "pointers_to_1"); + + if (!values_field || !indices_0_field || !indices_1_field || + !pointers_to_1_field) { + bsp_destroy_matrix_t(matrix); + return BSP_INVALID_STRUCTURE; + } + + bsp_error_t error = + matlab_to_bsp_array_allocator(values_field, &matrix->values, allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_matrix_t(matrix); + return error; + } + + error = matlab_to_bsp_array_allocator(indices_0_field, &matrix->indices_0, + allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_matrix_t(matrix); + return error; + } + + error = matlab_to_bsp_array_allocator(indices_1_field, &matrix->indices_1, + allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_matrix_t(matrix); + return error; + } + + error = matlab_to_bsp_array_allocator(pointers_to_1_field, + &matrix->pointers_to_1, allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_matrix_t(matrix); + return error; + } + + mxArray* nrows_field = mxGetField(mx_struct, 0, "nrows"); + mxArray* ncols_field = mxGetField(mx_struct, 0, "ncols"); + mxArray* nnz_field = mxGetField(mx_struct, 0, "nnz"); + mxArray* is_iso_field = mxGetField(mx_struct, 0, "is_iso"); + + if (!nrows_field || !ncols_field || !nnz_field || !is_iso_field) { + bsp_destroy_matrix_t(matrix); + return BSP_INVALID_STRUCTURE; + } + + matrix->nrows = (size_t) mxGetScalar(nrows_field); + matrix->ncols = (size_t) mxGetScalar(ncols_field); + matrix->nnz = (size_t) mxGetScalar(nnz_field); + matrix->is_iso = mxIsLogicalScalarTrue(is_iso_field); + + mxArray* format_field = mxGetField(mx_struct, 0, "format"); + if (!format_field || !mxIsChar(format_field)) { + bsp_destroy_matrix_t(matrix); + return BSP_INVALID_STRUCTURE; + } + + char* format_str = mxArrayToString(format_field); + if (!format_str) { + bsp_destroy_matrix_t(matrix); + return BSP_INVALID_STRUCTURE; + } + + matrix->format = bsp_get_matrix_format(format_str); + mxFree(format_str); + + if (matrix->format == BSP_INVALID_FORMAT) { + bsp_destroy_matrix_t(matrix); + return BSP_INVALID_FORMAT; + } + + mxArray* structure_field = mxGetField(mx_struct, 0, "structure"); + if (!structure_field || !mxIsChar(structure_field)) { + bsp_destroy_matrix_t(matrix); + return BSP_INVALID_STRUCTURE; + } + + char* structure_str = mxArrayToString(structure_field); + if (!structure_str) { + bsp_destroy_matrix_t(matrix); + return BSP_INVALID_STRUCTURE; + } + + matrix->structure = bsp_get_structure(structure_str); + mxFree(structure_str); + + if (matrix->structure == BSP_INVALID_STRUCTURE) { + matrix->structure = BSP_GENERAL; + } + + return BSP_SUCCESS; +} + +static inline bsp_error_t +bsp_matrix_copy_with_allocator(const bsp_matrix_t* input, bsp_matrix_t* output, + bsp_allocator_t allocator) { + bsp_construct_default_matrix_t_allocator(output, allocator); + + output->nrows = input->nrows; + output->ncols = input->ncols; + output->nnz = input->nnz; + output->is_iso = input->is_iso; + output->format = input->format; + output->structure = input->structure; + output->values.type = input->values.type; + output->indices_0.type = input->indices_0.type; + output->indices_1.type = input->indices_1.type; + output->pointers_to_1.type = input->pointers_to_1.type; + + if (input->values.size > 0) { + bsp_error_t error = bsp_construct_array_t_allocator( + &output->values, input->values.size, input->values.type, allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_matrix_t(output); + return error; + } + memcpy(output->values.data, input->values.data, + input->values.size * bsp_type_size(input->values.type)); + } + + if (input->indices_0.size > 0) { + bsp_error_t error = bsp_construct_array_t_allocator( + &output->indices_0, input->indices_0.size, input->indices_0.type, + allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_matrix_t(output); + return error; + } + memcpy(output->indices_0.data, input->indices_0.data, + input->indices_0.size * bsp_type_size(input->indices_0.type)); + } + + if (input->indices_1.size > 0) { + bsp_error_t error = bsp_construct_array_t_allocator( + &output->indices_1, input->indices_1.size, input->indices_1.type, + allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_matrix_t(output); + return error; + } + memcpy(output->indices_1.data, input->indices_1.data, + input->indices_1.size * bsp_type_size(input->indices_1.type)); + } + + if (input->pointers_to_1.size > 0) { + bsp_error_t error = bsp_construct_array_t_allocator( + &output->pointers_to_1, input->pointers_to_1.size, + input->pointers_to_1.type, allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_matrix_t(output); + return error; + } + memcpy(output->pointers_to_1.data, input->pointers_to_1.data, + input->pointers_to_1.size * + bsp_type_size(input->pointers_to_1.type)); + } + + return BSP_SUCCESS; +} + +static inline mxClassID get_mxClassID(bsp_type_t type) { + switch (type) { + case BSP_UINT8: + return mxUINT8_CLASS; + case BSP_UINT16: + return mxUINT16_CLASS; + case BSP_UINT32: + return mxUINT32_CLASS; + case BSP_UINT64: + return mxUINT64_CLASS; + case BSP_INT8: + return mxINT8_CLASS; + case BSP_INT16: + return mxINT16_CLASS; + case BSP_INT32: + return mxINT32_CLASS; + case BSP_INT64: + return mxINT64_CLASS; + case BSP_FLOAT32: + return mxSINGLE_CLASS; + case BSP_FLOAT64: + return mxDOUBLE_CLASS; + case BSP_BINT8: + return mxUINT8_CLASS; + case BSP_COMPLEX_FLOAT32: + return mxSINGLE_CLASS; + case BSP_COMPLEX_FLOAT64: + return mxDOUBLE_CLASS; + default: + return mxUNKNOWN_CLASS; + } +} + +static inline mxComplexity get_mxComplexity(bsp_type_t type) { + if (type == BSP_COMPLEX_FLOAT32 || type == BSP_COMPLEX_FLOAT64) { + return mxCOMPLEX; + } + return mxREAL; +} + +static inline mxArray* bsp_array_to_matlab(bsp_array_t* array) { + if (!array || array->data == NULL || array->size == 0) { + bsp_type_t type = array ? array->type : BSP_FLOAT64; + mxClassID class_id = get_mxClassID(type); + if (class_id == mxUNKNOWN_CLASS) { + class_id = mxDOUBLE_CLASS; + type = BSP_FLOAT64; + } + mxArray* empty_array = + mxCreateNumericMatrix(0, 1, class_id, get_mxComplexity(type)); + return empty_array; + } + + if (get_mxClassID(array->type) == mxUNKNOWN_CLASS) { + mexWarnMsgIdAndTxt("BinSparse:UnsupportedType", + "Unsupported array type %d, returning empty array", + (int) array->type); + mxArray* empty_array = mxCreateNumericMatrix(0, 1, mxDOUBLE_CLASS, mxREAL); + return empty_array; + } + + mxArray* mx_array = NULL; + + if ((array->allocator.malloc == bsp_matlab_allocator.malloc && + array->allocator.free == bsp_matlab_allocator.free) && + get_mxComplexity(array->type) == mxREAL) { + mx_array = mxCreateNumericMatrix(0, 1, get_mxClassID(array->type), + get_mxComplexity(array->type)); + + mxSetData(mx_array, array->data); + mxSetM(mx_array, array->size); + + array->data = NULL; + array->size = 0; + } else { + mx_array = mxCreateNumericMatrix(array->size, 1, get_mxClassID(array->type), + get_mxComplexity(array->type)); + + if (get_mxComplexity(array->type) == mxREAL) { + memcpy(mxGetData(mx_array), array->data, + array->size * bsp_type_size(array->type)); + } else { + if (array->type == BSP_COMPLEX_FLOAT32) { + float _Complex* in_data = (float _Complex*) array->data; + float* real_data = (float*) mxGetData(mx_array); + float* imag_data = (float*) mxGetImagData(mx_array); + for (size_t i = 0; i < array->size; i++) { + real_data[i] = crealf(in_data[i]); + imag_data[i] = cimagf(in_data[i]); + } + } else { + double _Complex* in_data = (double _Complex*) array->data; + double* real_data = mxGetPr(mx_array); + double* imag_data = mxGetPi(mx_array); + for (size_t i = 0; i < array->size; i++) { + real_data[i] = creal(in_data[i]); + imag_data[i] = cimag(in_data[i]); + } + } + } + } + + return mx_array; +} + +static inline mxArray* bsp_matrix_to_matlab_struct(bsp_matrix_t* matrix) { + const char* field_names[] = { + "values", "indices_0", "indices_1", "pointers_to_1", "nrows", + "ncols", "nnz", "is_iso", "format", "structure"}; + + mxArray* mx_struct = mxCreateStructMatrix(1, 1, 10, field_names); + + mxSetField(mx_struct, 0, "values", bsp_array_to_matlab(&matrix->values)); + mxSetField(mx_struct, 0, "indices_0", + bsp_array_to_matlab(&matrix->indices_0)); + mxSetField(mx_struct, 0, "indices_1", + bsp_array_to_matlab(&matrix->indices_1)); + mxSetField(mx_struct, 0, "pointers_to_1", + bsp_array_to_matlab(&matrix->pointers_to_1)); + + mxSetField(mx_struct, 0, "nrows", + mxCreateDoubleScalar((double) matrix->nrows)); + mxSetField(mx_struct, 0, "ncols", + mxCreateDoubleScalar((double) matrix->ncols)); + mxSetField(mx_struct, 0, "nnz", mxCreateDoubleScalar((double) matrix->nnz)); + mxSetField(mx_struct, 0, "is_iso", mxCreateLogicalScalar(matrix->is_iso)); + + mxSetField(mx_struct, 0, "format", + mxCreateString(bsp_get_matrix_format_string(matrix->format))); + mxSetField(mx_struct, 0, "structure", + mxCreateString(bsp_get_structure_string(matrix->structure))); + + return mx_struct; +} + +#endif diff --git a/bindings/matlab/test_binsparse_from_ssmc.m b/bindings/matlab/test_binsparse_from_ssmc.m new file mode 100644 index 0000000..3d86e20 --- /dev/null +++ b/bindings/matlab/test_binsparse_from_ssmc.m @@ -0,0 +1,33 @@ +function test_binsparse_from_ssmc() +% TEST_BINSPARSE_FROM_SSMC - Basic test for binsparse_from_ssmc MEX function + +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +fprintf('=== Testing binsparse_from_ssmc MEX function ===\n\n'); + +if exist('binsparse_from_ssmc', 'file') ~= 3 + error('binsparse_from_ssmc MEX function not found. Please compile it first.'); +end + +n = 4; +A = sparse([1 3], [2 4], [10 20], n, n); +Zeros = sparse([2 4], [1 3], [1 1], n, n); + +mat = binsparse_from_ssmc(A, Zeros, 'CSC'); + +assert(isstruct(mat)); +assert(mat.nrows == n && mat.ncols == n); +assert(strcmp(mat.format, 'CSC')); +assert(mat.nnz == nnz(A) + nnz(Zeros)); + +% Validate that explicit zeros were inserted +zero_values = mat.values(mat.values == 0); +if numel(zero_values) ~= nnz(Zeros) + error('Expected %d explicit zero values, got %d', nnz(Zeros), numel(zero_values)); +end + +fprintf('Test passed.\n'); + +end diff --git a/bindings/matlab/test_binsparse_minimize_roundtrip.m b/bindings/matlab/test_binsparse_minimize_roundtrip.m new file mode 100644 index 0000000..9c15e2b --- /dev/null +++ b/bindings/matlab/test_binsparse_minimize_roundtrip.m @@ -0,0 +1,67 @@ +function test_binsparse_minimize_roundtrip() +% TEST_BINSPARSE_MINIMIZE_ROUNDTRIP - Test SSMC conversion + type minimization + +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +fprintf('=== Testing binsparse_from_ssmc + binsparse_minimize_types ===\n\n'); + +if exist('binsparse_from_ssmc', 'file') ~= 3 + error('binsparse_from_ssmc MEX function not found. Please compile it first.'); +end + +if exist('binsparse_minimize_types', 'file') ~= 3 + error('binsparse_minimize_types MEX function not found. Please compile it first.'); +end + +if exist('binsparse_write', 'file') ~= 3 || exist('binsparse_read', 'file') ~= 3 + error('binsparse_write/binsparse_read MEX functions not found. Please compile them first.'); +end + +% Build a small matrix with explicit zeros pattern +n = 5; +A = sparse([1 3 5], [1 2 5], [1 2 3], n, n); +Zeros = sparse([2 4], [2 3], [1 1], n, n); + +mat = binsparse_from_ssmc(A, Zeros, 'CSC'); +mat_min = binsparse_minimize_types(mat); + +assert(strcmp(mat_min.format, 'CSC')); +assert(mat_min.nrows == n && mat_min.ncols == n); +assert(mat_min.nnz == mat.nnz); + +% Values should drop to single for exact float32 representable data +if ~strcmp(class(mat_min.values), 'single') + error('Expected values to minimize to single, got %s', class(mat_min.values)); +end + +% Index arrays should shrink to unsigned integer types +if ~strcmp(class(mat_min.indices_1), 'uint8') + error('Expected indices_1 to minimize to uint8, got %s', class(mat_min.indices_1)); +end + +if ~strcmp(class(mat_min.pointers_to_1), 'uint8') + error('Expected pointers_to_1 to minimize to uint8, got %s', class(mat_min.pointers_to_1)); +end + +% Roundtrip through binsparse_write/binsparse_read +out_file = [tempname() '.bsp.h5']; +cleanup_file = onCleanup(@() delete_if_exists(out_file)); + +binsparse_write(out_file, mat_min); +mat_read = binsparse_read(out_file); + +assert(mat_read.nrows == mat_min.nrows && mat_read.ncols == mat_min.ncols); +assert(mat_read.nnz == mat_min.nnz); +assert(strcmp(mat_read.format, mat_min.format)); + +fprintf('Test passed.\n'); + +end + +function delete_if_exists(filename) +if exist(filename, 'file') + delete(filename); +end +end diff --git a/bindings/matlab/test_binsparse_read.m b/bindings/matlab/test_binsparse_read.m new file mode 100644 index 0000000..bacd8d6 --- /dev/null +++ b/bindings/matlab/test_binsparse_read.m @@ -0,0 +1,70 @@ +function test_binsparse_read() +% TEST_BINSPARSE_READ - Test the binsparse_read MEX function +% +% This function demonstrates how to use the binsparse_read MEX function +% to read Binsparse format matrices into MATLAB/Octave. + +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +fprintf('=== Testing Binsparse Read Function ===\n\n'); + +% Check if the MEX function exists +if ~exist('binsparse_read', 'file') + error('binsparse_read MEX function not found. Build it first.'); +end + +try + % Test error handling with invalid inputs + fprintf('Test 1: Error handling\n'); + + try + matrix = binsparse_read(); % No arguments + fprintf('FAILED - Should have thrown error\n'); + catch ME + fprintf(' Correctly caught error for no arguments\n'); + end + + try + matrix = binsparse_read(123); % Non-string argument + fprintf('FAILED - Should have thrown error\n'); + catch ME + fprintf(' Correctly caught error for invalid filename type\n'); + end + + try + matrix = binsparse_read('nonexistent_file.h5'); % File doesn't exist + fprintf('FAILED - Should have thrown error\n'); + catch ME + fprintf(' Correctly caught error for nonexistent file\n'); + end + + fprintf(' Status: PASS\n\n'); + + fprintf('=== Basic Error Handling Tests Passed ===\n'); + fprintf('To test actual file reading, you need a valid Binsparse file.\n\n'); + + fprintf('Usage examples:\n'); + fprintf(' matrix = binsparse_read(''myfile.h5''); %% Read from HDF5\n'); + fprintf(' matrix = binsparse_read(''myfile.h5'', ''group''); %% Read from specific group\n\n'); + + fprintf('The returned matrix will have these fields:\n'); + fprintf(' matrix.values - Values array\n'); + fprintf(' matrix.indices_0 - First dimension indices\n'); + fprintf(' matrix.indices_1 - Second dimension indices\n'); + fprintf(' matrix.pointers_to_1 - Pointers for compressed formats\n'); + fprintf(' matrix.nrows - Number of rows\n'); + fprintf(' matrix.ncols - Number of columns\n'); + fprintf(' matrix.nnz - Number of non-zeros\n'); + fprintf(' matrix.is_iso - ISO matrix flag\n'); + fprintf(' matrix.format - Matrix format string\n'); + fprintf(' matrix.structure - Matrix structure string\n'); + +catch ME + fprintf('=== TEST FAILED ===\n'); + fprintf('Error: %s\n', ME.message); + rethrow(ME); +end + +end diff --git a/bindings/matlab/test_binsparse_roundtrip_dir.m b/bindings/matlab/test_binsparse_roundtrip_dir.m new file mode 100644 index 0000000..db96464 --- /dev/null +++ b/bindings/matlab/test_binsparse_roundtrip_dir.m @@ -0,0 +1,156 @@ +function test_binsparse_roundtrip_dir(root_dir, temp_dir) +% TEST_BINSPARSE_ROUNDTRIP_DIR - Round-trip binsparse files in a directory. +% +% This function scans a directory (recursively) for .h5 files, reads each +% with binsparse_read, writes to a temporary file with binsparse_write, then +% reads back and checks for equivalence. + +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +if nargin < 1 || isempty(root_dir) + error('Usage: test_binsparse_roundtrip_dir(root_dir, [temp_dir])'); +end + +if ~isfolder(root_dir) + error('Directory not found: %s', root_dir); +end + +if ~exist('binsparse_read', 'file') + error('binsparse_read MEX function not found. Build it first.'); +end + +if ~exist('binsparse_write', 'file') + error('binsparse_write MEX function not found. Build it first.'); +end + +if nargin < 2 || isempty(temp_dir) + temp_dir = ''; +end + +if ~isempty(temp_dir) && ~isfolder(temp_dir) + error('Temp directory not found: %s', temp_dir); +end + +files = list_h5_files(root_dir); +if isempty(files) + fprintf('No .h5 files found under %s\n', root_dir); + return; +end + +fprintf('Found %d .h5 files under %s\n', numel(files), root_dir); +failures = 0; + +for idx = 1:numel(files) + file_path = files{idx}; + fprintf('\n[%d/%d] %s\n', idx, numel(files), file_path); + + try + matrix = binsparse_read(file_path); + + if isempty(temp_dir) + temp_file = [tempname(), '.bsp.h5']; + else + temp_file = [tempname(temp_dir), '.bsp.h5']; + end + cleanup = onCleanup(@() cleanup_temp_file(temp_file)); + + binsparse_write(temp_file, matrix); + roundtrip = binsparse_read(temp_file); + + [ok, message] = compare_binsparse_structs(matrix, roundtrip); + if ok + fprintf(' OK\n'); + else + failures = failures + 1; + fprintf(' MISMATCH: %s\n', message); + end + + clear cleanup; + catch ME + failures = failures + 1; + fprintf(' ERROR: %s\n', ME.message); + end +end + +if failures == 0 + fprintf('\nAll %d files passed round-trip checks.\n', numel(files)); +else + fprintf('\n%d of %d files failed round-trip checks.\n', failures, numel(files)); +end + +end + +function files = list_h5_files(root_dir) +entries = dir(root_dir); +files = {}; + +for i = 1:numel(entries) + name = entries(i).name; + if entries(i).isdir + if strcmp(name, '.') || strcmp(name, '..') + continue; + end + sub_files = list_h5_files(fullfile(root_dir, name)); + if ~isempty(sub_files) + files = [files, sub_files]; %#ok + end + else + [~, ~, ext] = fileparts(name); + if strcmpi(ext, '.h5') + files{end + 1} = fullfile(root_dir, name); %#ok + end + end +end + +end + +function cleanup_temp_file(temp_file) +if exist(temp_file, 'file') + delete(temp_file); +end +end + +function [ok, message] = compare_binsparse_structs(a, b) +fields = {'values', 'indices_0', 'indices_1', 'pointers_to_1', ... + 'nrows', 'ncols', 'nnz', 'is_iso', 'format', 'structure'}; + +for i = 1:numel(fields) + field = fields{i}; + if ~isfield(a, field) + ok = false; + message = ['missing field in first matrix: ', field]; + return; + end + if ~isfield(b, field) + ok = false; + message = ['missing field in second matrix: ', field]; + return; + end + if ~compare_field(a.(field), b.(field)) + ok = false; + message = ['field mismatch: ', field]; + return; + end +end + +extra_a = setdiff(fieldnames(a), fields); +extra_b = setdiff(fieldnames(b), fields); +if ~isempty(extra_a) || ~isempty(extra_b) + ok = false; + message = 'field set mismatch'; + return; +end + +ok = true; +message = ''; +end + +function ok = compare_field(a, b) +if ischar(a) || isstring(a) || ischar(b) || isstring(b) + ok = strcmp(char(a), char(b)); +else + ok = isequaln(a, b); +end +end diff --git a/bindings/matlab/test_binsparse_write.m b/bindings/matlab/test_binsparse_write.m new file mode 100644 index 0000000..13ac3b2 --- /dev/null +++ b/bindings/matlab/test_binsparse_write.m @@ -0,0 +1,122 @@ +function test_binsparse_write() +% TEST_BINSPARSE_WRITE - Test the binsparse_write MEX function +% +% This function creates a simple test matrix and writes it to a temporary +% Binsparse file using the binsparse_write MEX function. + +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +fprintf('=== Testing binsparse_write MEX function ===\n\n'); + +% Check if binsparse_write exists +if ~exist('binsparse_write', 'file') + error('binsparse_write MEX function not found. Please compile it first.'); +end + +try + % Create a simple test matrix in the expected struct format + fprintf('Creating test matrix...\n'); + + % Simple 3x3 COO matrix: + % [1.0, 0.0, 2.0] + % [0.0, 3.0, 0.0] + % [4.0, 0.0, 5.0] + + matrix = struct(); + matrix.values = [1.0; 3.0; 2.0; 4.0; 5.0]; % Values in COO order + matrix.indices_0 = uint64([0; 1; 0; 2; 2]); % Row indices (0-based) + matrix.indices_1 = uint64([0; 1; 2; 0; 2]); % Column indices (0-based) + matrix.pointers_to_1 = uint64([]); % Empty for COO format + matrix.nrows = 3; + matrix.ncols = 3; + matrix.nnz = 5; + matrix.is_iso = false; + matrix.format = 'COO'; + matrix.structure = 'general'; + + fprintf('Matrix created:\n'); + fprintf(' Format: %s\n', matrix.format); + fprintf(' Size: %dx%d\n', matrix.nrows, matrix.ncols); + fprintf(' NNZ: %d\n', matrix.nnz); + fprintf(' Structure: %s\n', matrix.structure); + + % Create temporary filename + temp_file = tempname(); + temp_file = [temp_file, '.bsp.h5']; + + fprintf('\nWriting matrix to: %s\n', temp_file); + + % Test basic write + binsparse_write(temp_file, matrix); + fprintf('✓ Basic write successful\n'); + + % Test with group + binsparse_write(temp_file, matrix, 'test_group'); + fprintf('✓ Write with group successful\n'); + + % Test with group and JSON metadata + json_metadata = '{"test": "metadata", "created_by": "test_binsparse_write"}'; + binsparse_write(temp_file, matrix, 'test_group_json', json_metadata); + fprintf('✓ Write with group and JSON successful\n'); + + % Test with all parameters + binsparse_write(temp_file, matrix, 'test_group_full', json_metadata, 6); + fprintf('✓ Write with all parameters successful\n'); + + % Verify file exists and has reasonable size + if exist(temp_file, 'file') + info = dir(temp_file); + fprintf('✓ Output file created (size: %d bytes)\n', info.bytes); + else + error('Output file was not created'); + end + + % Test round-trip if binsparse_read is available + if exist('binsparse_read', 'file') + fprintf('\nTesting round-trip (write → read)...\n'); + + % Read back the matrix + read_matrix = binsparse_read(temp_file, 'test_group_full'); + + % Basic verification + if read_matrix.nrows == matrix.nrows && ... + read_matrix.ncols == matrix.ncols && ... + read_matrix.nnz == matrix.nnz + fprintf('✓ Round-trip dimensions match\n'); + else + warning('Round-trip dimensions mismatch'); + end + + if strcmp(read_matrix.format, matrix.format) && ... + strcmp(read_matrix.structure, matrix.structure) + fprintf('✓ Round-trip format/structure match\n'); + else + warning('Round-trip format/structure mismatch'); + end + else + fprintf('\nSkipping round-trip test (binsparse_read not available)\n'); + end + + % Clean up + if exist(temp_file, 'file') + delete(temp_file); + fprintf('✓ Temporary file cleaned up\n'); + end + + fprintf('\n=== All tests passed! ===\n'); + fprintf('binsparse_write is working correctly.\n\n'); + +catch ME + % Clean up on error + if exist('temp_file', 'var') && exist(temp_file, 'file') + delete(temp_file); + end + + fprintf('\n=== Test failed ===\n'); + fprintf('Error: %s\n', ME.message); + rethrow(ME); +end + +end diff --git a/bindings/matlab/test_bsp_matrix_struct.m b/bindings/matlab/test_bsp_matrix_struct.m new file mode 100644 index 0000000..a129593 --- /dev/null +++ b/bindings/matlab/test_bsp_matrix_struct.m @@ -0,0 +1,76 @@ +function test_bsp_matrix_struct() +% TEST_BSP_MATRIX_STRUCT - Test the Binsparse matrix struct functionality +% +% This function demonstrates and tests the basic MATLAB struct +% that mirrors the C bsp_matrix_t structure. + +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +fprintf('=== Testing Binsparse Matrix Struct ===\n\n'); + +try + % Test 1: Create empty matrix + fprintf('Test 1: Creating empty matrix\n'); + empty_matrix = bsp_matrix_create(); + fprintf('Empty matrix created successfully\n'); + bsp_matrix_info(empty_matrix); + fprintf('\n'); + + % Test 2: Create simple COO matrix + fprintf('Test 2: Creating simple COO matrix\n'); + % 3x3 identity matrix in COO format + values = [1.0, 1.0, 1.0]; + rows = uint64([0, 1, 2]); % 0-based indexing like C + cols = uint64([0, 1, 2]); % 0-based indexing like C + pointers = uint64([]); % Empty for COO format + + coo_matrix = bsp_matrix_create(values, rows, cols, pointers, ... + 3, 3, 3, false, 'COO', 'general'); + fprintf('COO matrix created successfully\n'); + bsp_matrix_info(coo_matrix); + fprintf('\n'); + + % Test 3: Create CSR matrix + fprintf('Test 3: Creating simple CSR matrix\n'); + % Same 3x3 identity in CSR format + csr_values = [1.0, 1.0, 1.0]; + csr_cols = uint64([0, 1, 2]); + csr_rows = uint64([]); % Not used in CSR + csr_ptrs = uint64([0, 1, 2, 3]); % Row pointers + + csr_matrix = bsp_matrix_create(csr_values, csr_rows, csr_cols, csr_ptrs, ... + 3, 3, 3, false, 'CSR', 'general'); + fprintf('CSR matrix created successfully\n'); + bsp_matrix_info(csr_matrix); + fprintf('\n'); + + % Test 4: Test field access + fprintf('Test 4: Testing field access\n'); + fprintf('Matrix format: %s\n', csr_matrix.format); + fprintf('Matrix structure: %s\n', csr_matrix.structure); + fprintf('Is ISO: %s\n', mat2str(csr_matrix.is_iso)); + fprintf('First value: %.1f\n', csr_matrix.values(1)); + fprintf('\n'); + + % Test 5: Test error handling + fprintf('Test 5: Testing error handling\n'); + try + invalid_matrix = bsp_matrix_create(1, 2, 3); % Wrong number of args + fprintf('FAILED - Should have thrown error\n'); + catch ME + fprintf('Successfully caught error: %s\n', ME.message); + end + fprintf('\n'); + + fprintf('=== All Tests Passed ===\n'); + fprintf('The Binsparse matrix struct is working correctly!\n'); + +catch ME + fprintf('=== TEST FAILED ===\n'); + fprintf('Error: %s\n', ME.message); + rethrow(ME); +end + +end diff --git a/bindings/matlab/test_convert_to_problem_struct.m b/bindings/matlab/test_convert_to_problem_struct.m new file mode 100644 index 0000000..46b8092 --- /dev/null +++ b/bindings/matlab/test_convert_to_problem_struct.m @@ -0,0 +1,161 @@ +function test_convert_to_problem_struct +%TEST_CONVERT_TO_PROBLEM_STRUCT test in-memory Binsparse Problem conversion + +% SPDX-FileCopyrightText: 2026 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +metadata = struct('name', 'Test/converter', 'title', 'converter test', ... + 'id', 7, 'date', '2026', 'author', 'Binsparse Developers', ... + 'ed', 'Binsparse Developers', 'kind', 'test', ... + 'notes', {{'line one'; 'line two'}}); + +rows = [0; 0; 1; 2; 2]; +cols = [0; 2; 1; 0; 2]; +values = [1; 0; 2; 3; 4]; +expected = sparse([1 2 3 3], [1 2 1 3], [1 2 3 4], 3, 3); +expected_zeros = sparse(1, 3, 1, 3, 3); + +formats = matrix_formats(rows, cols, values, 3, 3); +for k = 1:numel(formats) + raw = struct('metadata', metadata, 'A', formats{k}); + Problem = convert_to_problem_struct(raw); + assert(isequal(Problem.A, expected), ... + 'A mismatch for %s', formats{k}.format); + assert(isequal(Problem.Zeros, expected_zeros), ... + 'Zeros mismatch for %s', formats{k}.format); + assert(isequal(Problem.notes, char({'line one'; 'line two'}))); +end + +raw = struct('metadata', metadata, ... + 'A', make_matrix(0, rows, cols, [], 3, 3, 'COO', true, 'general')); +Problem = convert_to_problem_struct(raw); +assert(nnz(Problem.A) == 0); +assert(isequal(Problem.Zeros, sparse(rows + 1, cols + 1, 1, 3, 3))); + +lower = make_matrix([1; 2; 0; 3], [0; 1; 1; 2], [0; 0; 1; 2], ... + [], 3, 3, 'COO', false, 'symmetric_lower'); +raw = struct('metadata', metadata, 'A', lower); +Problem = convert_to_problem_struct(raw); +assert(isequal(Problem.A, sparse([1 2 1 3], [1 1 2 3], [1 2 2 3], 3, 3))); +assert(isequal(Problem.Zeros, sparse(2, 2, 1, 3, 3))); + +hermitian = make_matrix([1; 2+3i; 4], [0; 1; 1], [0; 0; 1], ... + [], 2, 2, 'COO', false, 'hermitian_lower'); +raw = struct('metadata', metadata, 'A', hermitian); +Problem = convert_to_problem_struct(raw); +assert(isequal(Problem.A, sparse([1 2-3i; 2+3i 4]))); + +skew = make_matrix([5; 0], [1; 1], [0; 1], ... + [], 2, 2, 'COO', false, 'skew_symmetric_lower'); +raw = struct('metadata', metadata, 'A', skew); +Problem = convert_to_problem_struct(raw); +assert(isequal(Problem.A, sparse([0 -5; 5 0]))); +assert(isequal(Problem.Zeros, sparse(2, 2, 1, 2, 2))); + +cvec = make_matrix([8; 0], [0; 2], [], [], ... + 3, 1, 'CVEC', false, 'general'); +raw = struct('metadata', metadata, 'A', cvec); +Problem = convert_to_problem_struct(raw); +assert(isequal(Problem.A, sparse(1, 1, 8, 3, 1))); +assert(isequal(Problem.Zeros, sparse(3, 1, 1, 3, 1))); + +raw = struct('metadata', metadata, ... + 'A', formats{1}, ... + 'b', dense_matrix([10; 20; 30], 3, 1, 'DVEC'), ... + 'x', dense_matrix([1; 3; 2; 4], 2, 2, 'DMATC')); +raw.aux.seq_1 = dense_matrix([5; 6], 2, 1, 'DVEC'); +raw.aux.seq_2 = dense_matrix([7; 8], 2, 1, 'DVEC'); +raw.aux.label = {'abc'; 'def'}; +Problem = convert_to_problem_struct(raw); +assert(isequal(Problem.b, [10; 20; 30])); +assert(isequal(Problem.x, [1 2; 3 4])); +assert(iscell(Problem.aux.seq) && numel(Problem.aux.seq) == 2); +assert(isequal(Problem.aux.seq{2}, [7; 8])); +assert(isequal(Problem.aux.label, char({'abc'; 'def'}))); + +dmat = dense_matrix([1; 2; 3; 4; 5; 6], 2, 3, 'DMAT'); +raw = struct('metadata', metadata, 'A', formats{1}, 'b', dmat); +Problem = convert_to_problem_struct(raw); +assert(isequal(Problem.b, [1 2 3; 4 5 6])); + +bad = formats{1}; +bad.indices_1([1 2]) = bad.indices_1([2 1]); +assert_throws(@() convert_to_problem_struct( ... + struct('metadata', metadata, 'A', bad)), 'BinSparse:InvalidMatrix'); + +bad = formats{3}; +bad.pointers_to_1 = uint64([0; 3; 2; 5]); +assert_throws(@() convert_to_problem_struct( ... + struct('metadata', metadata, 'A', bad)), 'BinSparse:InvalidMatrix'); + +bad = formats{5}; +bad.pointers_to_1 = uint64([0; 2; 2; 5]); +assert_throws(@() convert_to_problem_struct( ... + struct('metadata', metadata, 'A', bad)), 'BinSparse:InvalidMatrix'); + +bad = formats{1}; +bad.format = 'CUSTOM'; +assert_throws(@() convert_to_problem_struct( ... + struct('metadata', metadata, 'A', bad)), 'BinSparse:UnsupportedFormat'); + +bad = dense_matrix((1:6).', 2, 3, 'DMATC'); +assert_throws(@() convert_to_problem_struct( ... + struct('metadata', metadata, 'A', bad)), 'BinSparse:InvalidProblem'); + +bad = make_matrix(uint64(9007199254740992) + uint64(1), 0, 0, [], ... + 1, 1, 'COO', false, 'general'); +assert_throws(@() convert_to_problem_struct( ... + struct('metadata', metadata, 'A', bad)), 'BinSparse:InexactValue'); + +fprintf('test_convert_to_problem_struct: all tests passed\n'); + +end + +function assert_throws(f, identifier) +try + f(); +catch me + assert(strcmp(me.identifier, identifier), ... + 'Expected %s, got %s', identifier, me.identifier); + return; +end +error('Expected error %s was not thrown', identifier); +end + +function formats = matrix_formats(rows, cols, values, m, n) +formats = cell(7, 1); +formats{1} = make_matrix(values, rows, cols, [], m, n, 'COO', false, 'general'); +formats{2} = make_matrix(values([1 4 3 2 5]), cols([1 4 3 2 5]), ... + rows([1 4 3 2 5]), [], m, n, 'COOC', false, 'general'); +formats{3} = make_matrix(values, [], cols, [0; 2; 3; 5], ... + m, n, 'CSR', false, 'general'); +formats{4} = make_matrix(values([1 4 3 2 5]), [], rows([1 4 3 2 5]), ... + [0; 2; 3; 5], m, n, 'CSC', false, 'general'); +formats{5} = make_matrix(values, [0; 1; 2], cols, [0; 2; 3; 5], ... + m, n, 'DCSR', false, 'general'); +formats{6} = make_matrix(values([1 4 3 2 5]), [0; 1; 2], ... + rows([1 4 3 2 5]), [0; 2; 3; 5], m, n, 'DCSC', false, 'general'); +formats{7} = make_matrix(values, rows, cols, [], m, n, 'COOR', false, 'general'); +end + +function matrix = dense_matrix(values, m, n, format) +matrix = make_matrix(values, [], [], [], m, n, format, false, 'general'); +end + +function matrix = make_matrix(values, indices_0, indices_1, pointers, ... + m, n, format, is_iso, structure) +matrix = struct('values', values(:), ... + 'indices_0', uint64(indices_0(:)), ... + 'indices_1', uint64(indices_1(:)), ... + 'pointers_to_1', uint64(pointers(:)), ... + 'nrows', m, 'ncols', n, ... + 'nnz', max([numel(indices_0), numel(indices_1), numel(values)]), ... + 'is_iso', logical(is_iso), 'format', format, 'structure', structure); +if is_iso + matrix.nnz = numel(indices_0); +end +if any(strcmp(format, {'DVEC', 'DMAT', 'DMATC'})) + matrix.nnz = m * n; +end +end diff --git a/bindings/matlab/test_generate_bsp_from_ssmc.m b/bindings/matlab/test_generate_bsp_from_ssmc.m new file mode 100644 index 0000000..9771b24 --- /dev/null +++ b/bindings/matlab/test_generate_bsp_from_ssmc.m @@ -0,0 +1,182 @@ +function test_generate_bsp_from_ssmc() +% TEST_GENERATE_BSP_FROM_SSMC - End-to-end test for generate_bsp_from_ssmc + +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +fprintf('=== Testing generate_bsp_from_ssmc ===\n\n'); + +required = {'binsparse_from_ssmc', 'binsparse_minimize_types', ... + 'binsparse_write', 'binsparse_read', ... + 'binsparse_write_string_dataset', 'generate_bsp_from_ssmc'}; +for i = 1:numel(required) + if exist(required{i}, 'file') ~= 3 && exist(required{i}, 'file') ~= 2 + error('%s not found. Please compile MEX functions and ensure the .m file is on path.', required{i}); + end +end + +% Build synthetic problem +Problem = struct(); +Problem.name = 'Test/Small'; +Problem.title = 'Small test'; +Problem.id = 7; +Problem.date = '2026-06-03'; +Problem.author = 'Binsparse Developers'; +Problem.ed = 'Binsparse Developers'; +Problem.kind = 'test matrix'; +Problem.notes = ['first note '; 'second note ']; +Problem.A = sparse([1 3 4], [1 2 4], [5 6 7], 4, 4); +Problem.Zeros = sparse([2 4], [2 1], [1 1], 4, 4); +Problem.b = [10; 20; 30; 40]; +Problem.x = [1 2; 3 4; 5 6; 7 8]; +Problem.aux = struct(); +Problem.aux.c = [1; 2; 3]; +Problem.aux.D = [1 0 2; 3 4 5]; +Problem.aux.S = sparse([1 2], [2 3], [9 8], 3, 3); +Problem.aux.note = ['hello '; 'there ']; +Problem.aux.tags = {'alpha'; 'beta'}; + +problem = struct('Problem', Problem); + +out_file = [tempname() '.bsp.h5']; +cleanup_file = onCleanup(@() delete_if_exists(out_file)); + +format = 'COO'; +compression_level = 0; + +generate_bsp_from_ssmc(problem, out_file, format, compression_level); + +% Read primary +primary_bsp = binsparse_read(out_file); +primary_mat = bsp_to_matlab(primary_bsp); +expected_primary = full(Problem.A); +assert(matrices_equal(primary_mat, expected_primary), 'Primary matrix mismatch'); + +% Read aux and x/b +check_dense_group(out_file, 'b', Problem.b(:)); +check_dense_group(out_file, 'x', Problem.x); +check_dense_group(out_file, 'c', Problem.aux.c(:)); +check_dense_group(out_file, 'D', Problem.aux.D); + +aux_sparse = binsparse_read(out_file, 'S'); +aux_sparse_mat = bsp_to_matlab(aux_sparse); +expected_sparse = full(Problem.aux.S); +assert(matrices_equal(aux_sparse_mat, expected_sparse), 'Aux sparse matrix mismatch'); + +check_string_dataset(out_file, 'note', {'hello '; 'there '}); +check_string_dataset(out_file, 'tags', {'alpha'; 'beta'}); + +json = h5readatt(out_file, '/', 'binsparse'); +assert(contains(json, '"metadata"'), 'Primary metadata not nested'); +assert(~isempty(regexp(json, '"id"\s*:\s*7', 'once')), ... + 'Primary metadata id is not a JSON number'); +assert(~contains(json, '"ssmc_metadata"'), 'Unexpected legacy metadata key'); + +fprintf('Test passed.\n'); + +end + +function check_dense_group(filename, group, expected) + bsp = binsparse_read(filename, group); + actual = bsp_to_matlab(bsp); + assert(matrices_equal(actual, expected), ... + 'Group "%s" mismatch', group); +end + +function check_string_dataset(filename, name, expected) + actual = h5read(filename, ['/' name]); + actual = as_cellstr(actual); + assert(isequal(actual, expected), 'String dataset "%s" mismatch', name); +end + +function value = as_cellstr(value) + if iscell(value) + value = value(:); + elseif isstring(value) + value = cellstr(value(:)); + elseif ischar(value) + value = cellstr(value); + else + error('Unexpected string dataset value type'); + end +end + +function ok = matrices_equal(a, b) + if ~isequal(size(a), size(b)) + ok = false; + return; + end + a = double(a); + b = double(b); + diff = a - b; + ok = all(diff(:) == 0); +end + +function mat = bsp_to_matlab(bsp) + fmt = upper(bsp.format); + switch fmt + case 'COO' + rows = double(bsp.indices_0) + 1; + cols = double(bsp.indices_1) + 1; + vals = bsp.values; + mat = sparse(rows, cols, vals, bsp.nrows, bsp.ncols); + mat = full(mat); + case 'CSC' + colptr = double(bsp.pointers_to_1); + rowind = double(bsp.indices_1); + vals = bsp.values; + ncols = bsp.ncols; + nrows = bsp.nrows; + rows = []; + cols = []; + v = []; + for j = 1:ncols + start_idx = colptr(j) + 1; + end_idx = colptr(j + 1); + if end_idx >= start_idx + idx = start_idx:end_idx; + rows = [rows; rowind(idx) + 1]; + cols = [cols; j * ones(numel(idx), 1)]; + v = [v; vals(idx)]; + end + end + mat = sparse(rows, cols, v, nrows, ncols); + mat = full(mat); + case 'CSR' + rowptr = double(bsp.pointers_to_1); + colind = double(bsp.indices_1); + vals = bsp.values; + nrows = bsp.nrows; + ncols = bsp.ncols; + rows = []; + cols = []; + v = []; + for i = 1:nrows + start_idx = rowptr(i) + 1; + end_idx = rowptr(i + 1); + if end_idx >= start_idx + idx = start_idx:end_idx; + rows = [rows; i * ones(numel(idx), 1)]; + cols = [cols; colind(idx) + 1]; + v = [v; vals(idx)]; + end + end + mat = sparse(rows, cols, v, nrows, ncols); + mat = full(mat); + case 'DMAT' + mat = reshape(bsp.values, [bsp.ncols, bsp.nrows]).'; + case 'DMATC' + mat = reshape(bsp.values, [bsp.nrows, bsp.ncols]); + case 'DVEC' + mat = reshape(bsp.values, [bsp.nrows, 1]); + otherwise + error('Unsupported format in test: %s', fmt); + end +end + +function delete_if_exists(filename) + if exist(filename, 'file') + delete(filename); + end +end diff --git a/bindings/matlab/test_write_binsparse_from_matlab.m b/bindings/matlab/test_write_binsparse_from_matlab.m new file mode 100644 index 0000000..a663910 --- /dev/null +++ b/bindings/matlab/test_write_binsparse_from_matlab.m @@ -0,0 +1,102 @@ +function test_write_binsparse_from_matlab() +% TEST_WRITE_BINSPARSE_FROM_MATLAB - End-to-end test for the SSMC writer MEX + +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +fprintf('=== Testing write_binsparse_from_matlab MEX function ===\n\n'); + +required = {'write_binsparse_from_matlab', 'binsparse_read'}; +for i = 1:numel(required) + if exist(required{i}, 'file') ~= 3 + error('%s MEX function not found. Please compile it first.', required{i}); + end +end + +try + write_binsparse_from_matlab(); + error('Expected insufficient-argument error'); +catch ME + assert(strcmp(ME.identifier, 'BinSparse:InvalidArgs'), ... + 'Unexpected identifier: %s', ME.identifier); +end + +try + write_binsparse_from_matlab(42, 'invalid.bsp.h5'); + error('Expected invalid-problem error'); +catch ME + assert(strcmp(ME.identifier, 'BinSparse:InvalidProblem'), ... + 'Unexpected identifier: %s', ME.identifier); +end + +Problem = struct(); +Problem.name = 'test_basic_matrix'; +Problem.title = 'Basic Test Matrix'; +Problem.kind = 'test matrix'; +Problem.A = sparse([1 3], [1 2], [5 7], 3, 3); +Problem.Zeros = sparse(2, 3, 1, 3, 3); +Problem.b = [1; 2; 3]; +Problem.aux = struct(); +Problem.aux.D = [10 20; 30 40]; + +out_file = [tempname() '.bsp.h5']; +cleanup_file = onCleanup(@() delete_if_exists(out_file)); + +write_binsparse_from_matlab(struct('Problem', Problem), out_file, 'COO', [], 0); + +primary = bsp_to_matlab(binsparse_read(out_file)); +assert(isequal(double(primary), full(Problem.A)), 'Primary matrix mismatch'); + +b = bsp_to_matlab(binsparse_read(out_file, 'b')); +assert(isequal(double(b), Problem.b), 'b vector mismatch'); + +D = bsp_to_matlab(binsparse_read(out_file, 'D')); +assert(isequal(double(D), Problem.aux.D), 'aux.D mismatch'); + +legacy_out_file = [tempname() '.bsp.h5']; +cleanup_legacy = onCleanup(@() delete_if_exists(legacy_out_file)); +write_binsparse_from_matlab(Problem, legacy_out_file, 'legacy_group', '{}', 0); +assert(exist(legacy_out_file, 'file') == 2, 'Legacy call did not create file'); + +fprintf('Test passed.\n'); + +end + +function mat = bsp_to_matlab(bsp) + fmt = upper(bsp.format); + switch fmt + case 'COO' + rows = double(bsp.indices_0) + 1; + cols = double(bsp.indices_1) + 1; + vals = bsp.values; + mat = full(sparse(rows, cols, vals, bsp.nrows, bsp.ncols)); + case 'CSC' + colptr = double(bsp.pointers_to_1); + rowind = double(bsp.indices_1); + rows = []; + cols = []; + vals = []; + for j = 1:bsp.ncols + idx = (colptr(j) + 1):colptr(j + 1); + rows = [rows; rowind(idx) + 1]; + cols = [cols; j * ones(numel(idx), 1)]; + vals = [vals; bsp.values(idx)]; + end + mat = full(sparse(rows, cols, vals, bsp.nrows, bsp.ncols)); + case 'DMAT' + mat = reshape(bsp.values, [bsp.ncols, bsp.nrows]).'; + case 'DMATC' + mat = reshape(bsp.values, [bsp.nrows, bsp.ncols]); + case 'DVEC' + mat = reshape(bsp.values, [bsp.nrows, 1]); + otherwise + error('Unsupported format in test: %s', fmt); + end +end + +function delete_if_exists(filename) +if exist(filename, 'file') + delete(filename); +end +end diff --git a/bindings/matlab/write_binsparse_from_matlab.c b/bindings/matlab/write_binsparse_from_matlab.c new file mode 100644 index 0000000..69951b8 --- /dev/null +++ b/bindings/matlab/write_binsparse_from_matlab.c @@ -0,0 +1,116 @@ +/* + * SPDX-FileCopyrightText: 2024 Binsparse Developers + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +/** + * write_binsparse_from_matlab.c - Write a SuiteSparse Matrix Collection + * problem struct to Binsparse format. + * + * This MEX entry point intentionally delegates SSMC-specific handling to the + * MATLAB implementation in generate_bsp_from_ssmc.m. That path writes the + * primary matrix, explicit zeros, metadata, aux matrices, b, and x. + * + * Usage in MATLAB/Octave: + * write_binsparse_from_matlab(problem_struct, filename) + * write_binsparse_from_matlab(problem_struct, filename, format) + * write_binsparse_from_matlab(problem_struct, filename, format, [], + * compression) + * + * For compatibility with the older scaffold, a non-format third string is + * accepted and ignored as an HDF5 group name; generate_bsp_from_ssmc writes + * SSMC files at the root with aux entries in named groups. + */ + +#include "mex.h" +#include +#include + +static bool is_supported_sparse_format(const char* text) { + return strcasecmp(text, "CSC") == 0 || strcasecmp(text, "CSR") == 0 || + strcasecmp(text, "COO") == 0 || strcasecmp(text, "COOR") == 0; +} + +static mxArray* make_default_format(void) { + return mxCreateString("COO"); +} + +void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { + (void) plhs; + + if (nrhs < 2 || nrhs > 5) { + mexErrMsgIdAndTxt( + "BinSparse:InvalidArgs", + "Usage: write_binsparse_from_matlab(problem_struct, filename [, format " + "[, json_metadata [, compression_level]]])"); + } + + if (nlhs > 0) { + mexErrMsgIdAndTxt("BinSparse:TooManyOutputs", + "No output arguments expected"); + } + + if (!mxIsStruct(prhs[0])) { + mexErrMsgIdAndTxt("BinSparse:InvalidProblem", + "First argument must be a struct"); + } + + if (!mxIsChar(prhs[1])) { + mexErrMsgIdAndTxt("BinSparse:InvalidFilename", + "Filename must be a character vector"); + } + + mxArray* format = NULL; + if (nrhs >= 3 && !mxIsEmpty(prhs[2])) { + if (!mxIsChar(prhs[2])) { + mexErrMsgIdAndTxt("BinSparse:InvalidFormat", + "Format must be a character vector"); + } + + char* format_text = mxArrayToString(prhs[2]); + if (!format_text) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to read format string"); + } + + if (is_supported_sparse_format(format_text)) { + format = mxCreateString(format_text); + } + mxFree(format_text); + } + + if (!format) { + format = make_default_format(); + } + + mxArray* compression = NULL; + if (nrhs >= 5 && !mxIsEmpty(prhs[4])) { + if (!mxIsNumeric(prhs[4]) || mxIsComplex(prhs[4]) || + mxGetNumberOfElements(prhs[4]) != 1) { + mxDestroyArray(format); + mexErrMsgIdAndTxt("BinSparse:InvalidCompression", + "Compression level must be a real numeric scalar"); + } + compression = mxDuplicateArray(prhs[4]); + } else if (nrhs >= 4 && mxIsNumeric(prhs[3]) && !mxIsEmpty(prhs[3])) { + compression = mxDuplicateArray(prhs[3]); + } else { + compression = mxCreateDoubleScalar(0.0); + } + + mxArray* rhs[4]; + rhs[0] = (mxArray*) prhs[0]; + rhs[1] = (mxArray*) prhs[1]; + rhs[2] = format; + rhs[3] = compression; + + int status = mexCallMATLAB(0, NULL, 4, rhs, "generate_bsp_from_ssmc"); + + mxDestroyArray(format); + mxDestroyArray(compression); + + if (status != 0) { + mexErrMsgIdAndTxt("BinSparse:WriteFailed", "generate_bsp_from_ssmc failed"); + } +} diff --git a/examples/benchmark_read.c b/examples/benchmark_read.c index c1ef342..0a1894a 100644 --- a/examples/benchmark_read.c +++ b/examples/benchmark_read.c @@ -7,6 +7,7 @@ #include #include #include +#include double gettime() { struct timespec time; @@ -79,8 +80,9 @@ int main(int argc, char** argv) { // If running warm cache experiments, read once to warm cache. if (!cold_cache) { - bsp_matrix_t mat = bsp_read_matrix(file_name, NULL); - bsp_destroy_matrix_t(mat); + bsp_matrix_t mat; + BSP_CHECK(bsp_read_matrix(&mat, file_name, NULL)); + bsp_destroy_matrix_t(&mat); } for (size_t i = 0; i < num_trials; i++) { @@ -88,11 +90,12 @@ int main(int argc, char** argv) { flush_cache(); } double begin = gettime(); - bsp_matrix_t mat = bsp_read_matrix(file_name, NULL); + bsp_matrix_t mat; + BSP_CHECK(bsp_read_matrix(&mat, file_name, NULL)); double end = gettime(); durations[i] = end - begin; nbytes = bsp_matrix_nbytes(mat); - bsp_destroy_matrix_t(mat); + bsp_destroy_matrix_t(&mat); double gbytes = ((double) nbytes) / 1024 / 1024 / 1024; double gbytes_s = gbytes / durations[i]; diff --git a/examples/benchmark_read_parallel.c b/examples/benchmark_read_parallel.c index b485938..24a5597 100644 --- a/examples/benchmark_read_parallel.c +++ b/examples/benchmark_read_parallel.c @@ -7,6 +7,7 @@ #include #include #include +#include double gettime() { struct timespec time; @@ -82,8 +83,9 @@ int main(int argc, char** argv) { // If running warm cache experiments, read once to warm cache. if (!cold_cache) { - bsp_matrix_t mat = bsp_read_matrix_parallel(file_name, NULL, num_threads); - bsp_destroy_matrix_t(mat); + bsp_matrix_t mat; + BSP_CHECK(bsp_read_matrix_parallel(&mat, file_name, NULL, num_threads)); + bsp_destroy_matrix_t(&mat); } for (size_t i = 0; i < num_trials; i++) { @@ -92,12 +94,13 @@ int main(int argc, char** argv) { } fflush(stdout); double begin = gettime(); - bsp_matrix_t mat = bsp_read_matrix_parallel(file_name, NULL, num_threads); + bsp_matrix_t mat; + BSP_CHECK(bsp_read_matrix_parallel(&mat, file_name, NULL, num_threads)); double end = gettime(); durations[i] = end - begin; nbytes = bsp_matrix_nbytes(mat); - bsp_destroy_matrix_t(mat); + bsp_destroy_matrix_t(&mat); double gbytes = ((double) nbytes) / 1024 / 1024 / 1024; double gbytes_s = gbytes / durations[i]; diff --git a/examples/benchmark_write.c b/examples/benchmark_write.c index 8af5881..eaed5ab 100644 --- a/examples/benchmark_write.c +++ b/examples/benchmark_write.c @@ -5,8 +5,10 @@ */ #include +#include #include #include +#include double gettime() { struct timespec time; @@ -69,9 +71,10 @@ void flush_writes() { } void delete_file(const char* file_name) { - char command[2048]; - snprintf(command, 2047, "rm %s", file_name); - int rv = system(command); + int rv = remove(file_name); + if (rv != 0) { + perror("delete_file"); + } } int main(int argc, char** argv) { @@ -97,7 +100,8 @@ int main(int argc, char** argv) { double durations[num_trials]; - bsp_matrix_t mat = bsp_read_matrix(file_name, NULL); + bsp_matrix_t mat; + BSP_CHECK(bsp_read_matrix(&mat, file_name, NULL)); size_t nbytes = bsp_matrix_nbytes(mat); char output_filename[2048]; @@ -124,7 +128,8 @@ int main(int argc, char** argv) { printf("Writing to file %s\n", output_filename); double begin = gettime(); - bsp_write_matrix(output_filename, mat, NULL, NULL, compression_level); + BSP_CHECK( + bsp_write_matrix(output_filename, mat, NULL, NULL, compression_level)); if (flush_each_write) { flush_writes(); diff --git a/examples/bsp-ls.c b/examples/bsp-ls.c index 8ecdd94..7fb0870 100644 --- a/examples/bsp-ls.c +++ b/examples/bsp-ls.c @@ -4,7 +4,7 @@ * SPDX-License-Identifier: BSD-3-Clause */ -#include +#include herr_t visit_group(hid_t loc_id, const char* name, const H5L_info_t* linfo, void* opdata); @@ -18,7 +18,13 @@ void print_group_info(hid_t g, const char* name) { H5E_END_TRY; if (bsp_json != H5I_INVALID_HID) { - char* json_string = bsp_read_attribute(g, "binsparse"); + char* json_string; + bsp_error_t error = bsp_read_attribute(&json_string, g, "binsparse"); + if (error != BSP_SUCCESS) { + printf("Error reading binsparse attribute: %s\n", + bsp_get_error_string(error)); + return; + } cJSON* j = cJSON_Parse(json_string); diff --git a/examples/bsp2mtx.c b/examples/bsp2mtx.c index 6e3f483..5a8f349 100644 --- a/examples/bsp2mtx.c +++ b/examples/bsp2mtx.c @@ -10,8 +10,7 @@ int main(int argc, char** argv) { if (argc < 3) { - printf( - "usage: ./bsp2mtx [inputfile_name.mtx] [outputfile_name.bsp.hdf5]\n"); + printf("usage: ./bsp2mtx [inputfile_name.bsp.h5] [outputfile_name.mtx]\n"); return 1; } @@ -19,8 +18,9 @@ int main(int argc, char** argv) { char* output_fname = argv[2]; printf(" === Reading file... ===\n"); - bsp_matrix_t matrix = bsp_read_matrix(input_fname, NULL); - printf(" === Done writing. ===\n"); + bsp_matrix_t matrix; + BSP_CHECK(bsp_read_matrix(&matrix, input_fname, NULL)); + printf(" === Done reading. ===\n"); if (matrix.format != BSP_COO) { matrix = bsp_convert_matrix(matrix, BSP_COO); } @@ -29,7 +29,7 @@ int main(int argc, char** argv) { bsp_mmwrite(output_fname, matrix); printf(" === Done writing. ===\n"); - bsp_destroy_matrix_t(matrix); + bsp_destroy_matrix_t(&matrix); return 0; } diff --git a/examples/check_equivalence.c b/examples/check_equivalence.c index 0acea42..641ec8c 100644 --- a/examples/check_equivalence.c +++ b/examples/check_equivalence.c @@ -110,8 +110,9 @@ int main(int argc, char** argv) { printf("Matrix 2: %s and %s\n", info2.fname, (info2.dataset == NULL) ? "root" : info2.dataset); - bsp_matrix_t matrix1 = bsp_read_matrix(info1.fname, info1.dataset); - bsp_matrix_t matrix2 = bsp_read_matrix(info2.fname, info2.dataset); + bsp_matrix_t matrix1, matrix2; + BSP_CHECK(bsp_read_matrix(&matrix1, info1.fname, info1.dataset)); + BSP_CHECK(bsp_read_matrix(&matrix2, info2.fname, info2.dataset)); bool perform_suitesparse_declamping = true; if (perform_suitesparse_declamping && @@ -128,13 +129,13 @@ int main(int argc, char** argv) { if (matrix1.format != matrix2.format) { if (matrix1.format != BSP_COOR) { bsp_matrix_t intermediate = bsp_convert_matrix(matrix1, BSP_COOR); - bsp_destroy_matrix_t(matrix1); + bsp_destroy_matrix_t(&matrix1); matrix1 = intermediate; } if (matrix2.format != BSP_COOR) { bsp_matrix_t intermediate = bsp_convert_matrix(matrix2, BSP_COOR); - bsp_destroy_matrix_t(matrix2); + bsp_destroy_matrix_t(&matrix2); matrix2 = intermediate; } } @@ -202,6 +203,11 @@ int main(int argc, char** argv) { return 9; } + bsp_destroy_matrix_t(&matrix1); + bsp_destroy_matrix_t(&matrix2); + bsp_destroy_fdataset_info_t(&info1); + bsp_destroy_fdataset_info_t(&info2); + printf("The files are equivalent.\n"); printf("OK!\n"); diff --git a/examples/check_equivalence_parallel.c b/examples/check_equivalence_parallel.c index b773cf8..ecdcaf0 100644 --- a/examples/check_equivalence_parallel.c +++ b/examples/check_equivalence_parallel.c @@ -117,9 +117,10 @@ int main(int argc, char** argv) { fflush(stdout); - bsp_matrix_t matrix1 = bsp_read_matrix(info1.fname, info1.dataset); - bsp_matrix_t matrix2 = - bsp_read_matrix_parallel(info2.fname, info2.dataset, num_threads); + bsp_matrix_t matrix1, matrix2; + BSP_CHECK(bsp_read_matrix(&matrix1, info1.fname, info1.dataset)); + BSP_CHECK(bsp_read_matrix_parallel(&matrix2, info2.fname, info2.dataset, + num_threads)); bool perform_suitesparse_declamping = true; if (perform_suitesparse_declamping && @@ -136,13 +137,13 @@ int main(int argc, char** argv) { if (matrix1.format != matrix2.format) { if (matrix1.format != BSP_COOR) { bsp_matrix_t intermediate = bsp_convert_matrix(matrix1, BSP_COOR); - bsp_destroy_matrix_t(matrix1); + bsp_destroy_matrix_t(&matrix1); matrix1 = intermediate; } if (matrix2.format != BSP_COOR) { bsp_matrix_t intermediate = bsp_convert_matrix(matrix2, BSP_COOR); - bsp_destroy_matrix_t(matrix2); + bsp_destroy_matrix_t(&matrix2); matrix2 = intermediate; } } diff --git a/examples/mtx2bsp.c b/examples/mtx2bsp.c index 07db200..c6f8083 100644 --- a/examples/mtx2bsp.c +++ b/examples/mtx2bsp.c @@ -4,7 +4,7 @@ * SPDX-License-Identifier: BSD-3-Clause */ -#include +#include #include #include @@ -99,11 +99,25 @@ int main(int argc, char** argv) { } bsp_mm_metadata m = bsp_mmread_metadata(input_fname); + if (!bsp_mm_metadata_is_valid(m)) { + fprintf(stderr, + "error: unable to read Matrix Market metadata from \"%s\".\n", + input_fname); + bsp_destroy_mm_metadata(&m); + bsp_destroy_fdataset_info_t(&info2); + return 1; + } bsp_matrix_format_t format = BSP_COOR; if (format_name != NULL) { format = bsp_get_matrix_format(format_name); - assert(format != 0); + if (format == BSP_INVALID_FORMAT) { + fprintf(stderr, "error: unsupported Binsparse format \"%s\".\n", + format_name); + bsp_destroy_mm_metadata(&m); + bsp_destroy_fdataset_info_t(&info2); + return 1; + } } printf("%lu x %lu matrix with %lu nonzeros.\n", m.nrows, m.ncols, m.nnz); @@ -151,7 +165,7 @@ int main(int argc, char** argv) { if (format != BSP_COOR) { begin = gettime(); bsp_matrix_t converted_matrix = bsp_convert_matrix(matrix, format); - bsp_destroy_matrix_t(matrix); + bsp_destroy_matrix_t(&matrix); matrix = converted_matrix; end = gettime(); duration = end - begin; @@ -163,14 +177,17 @@ int main(int argc, char** argv) { printf(" === Writing to %s... ===\n", output_fname); begin = gettime(); - bsp_write_matrix(output_fname, matrix, group_name, user_json, - compression_level); + BSP_CHECK(bsp_write_matrix_cjson(output_fname, matrix, group_name, user_json, + compression_level)); end = gettime(); duration = end - begin; printf("%lf seconds writing Binsparse file...\n", duration); printf(" === Done writing. ===\n"); - bsp_destroy_matrix_t(matrix); + bsp_destroy_matrix_t(&matrix); + bsp_destroy_mm_metadata(&m); + bsp_destroy_fdataset_info_t(&info2); + cJSON_Delete(user_json); return 0; } diff --git a/examples/simple_matrix_read.c b/examples/simple_matrix_read.c index acbdabc..8ca1411 100644 --- a/examples/simple_matrix_read.c +++ b/examples/simple_matrix_read.c @@ -9,7 +9,8 @@ int main(int argc, char** argv) { const char* file_name = "test.hdf5"; - bsp_matrix_t mat = bsp_read_matrix(file_name, NULL); + bsp_matrix_t mat; + BSP_CHECK(bsp_read_matrix(&mat, file_name, NULL)); if (mat.format == BSP_COO) { float* values = (float*) mat.values.data; diff --git a/examples/simple_matrix_write.c b/examples/simple_matrix_write.c index 5ee6f3f..852c7e1 100644 --- a/examples/simple_matrix_write.c +++ b/examples/simple_matrix_write.c @@ -21,7 +21,7 @@ int main(int argc, char** argv) { printf("%d, %d: %f\n", rowind[i], colind[i], values[i]); } - bsp_write_matrix("test.hdf5", mat, NULL, NULL, 9); + BSP_CHECK(bsp_write_matrix("test.hdf5", mat, NULL, NULL, 9)); return 0; } diff --git a/examples/simple_read.c b/examples/simple_read.c index a4ba9bb..18d2eac 100644 --- a/examples/simple_read.c +++ b/examples/simple_read.c @@ -4,14 +4,15 @@ * SPDX-License-Identifier: BSD-3-Clause */ -#include +#include int main(int argc, char** argv) { char* file_name = (char*) "test.hdf5"; hid_t f = H5Fopen(file_name, H5F_ACC_RDWR, H5P_DEFAULT); - bsp_array_t array = bsp_read_array(f, "test"); + bsp_array_t array; + BSP_CHECK(bsp_read_array(&array, f, "test")); int* values = (int*) array.data; @@ -19,7 +20,7 @@ int main(int argc, char** argv) { printf("%lu: %d\n", i, values[i]); } - bsp_destroy_array_t(array); + bsp_destroy_array_t(&array); H5Fclose(f); return 0; diff --git a/examples/simple_write.c b/examples/simple_write.c index deb47d8..0a51b49 100644 --- a/examples/simple_write.c +++ b/examples/simple_write.c @@ -4,14 +4,15 @@ * SPDX-License-Identifier: BSD-3-Clause */ -#include +#include int main(int argc, char** argv) { const char* file_name = "test.hdf5"; hid_t f = H5Fcreate(file_name, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); - bsp_array_t array = bsp_construct_array_t(1000, BSP_INT32); + bsp_array_t array; + bsp_construct_array_t(&array, 1000, BSP_INT32); int* values = (int*) array.data; @@ -22,7 +23,7 @@ int main(int argc, char** argv) { bsp_write_array(f, "test", array, 0); H5Fclose(f); - bsp_destroy_array_t(array); + bsp_destroy_array_t(&array); return 0; } diff --git a/include/binsparse/array.h b/include/binsparse/array.h index 8ac69af..0570d10 100644 --- a/include/binsparse/array.h +++ b/include/binsparse/array.h @@ -6,9 +6,9 @@ #pragma once -#include #include #include +#include #include #include #include @@ -20,71 +20,91 @@ typedef struct bsp_array_t { bsp_allocator_t allocator; } bsp_array_t; -static inline bsp_array_t bsp_construct_default_array_t() { - bsp_array_t array; - array.data = NULL; - array.size = 0; - array.allocator = bsp_default_allocator; - return array; +static inline bsp_error_t +bsp_construct_default_array_t_allocator(bsp_array_t* array, + bsp_allocator_t allocator) { + array->data = NULL; + array->size = 0; + array->type = BSP_INVALID_TYPE; + array->allocator = allocator; + return BSP_SUCCESS; } -static inline bsp_array_t bsp_construct_array_t(size_t size, bsp_type_t type) { +static inline bsp_error_t bsp_construct_default_array_t(bsp_array_t* array) { + return bsp_construct_default_array_t_allocator(array, bsp_default_allocator); +} + +static inline bsp_error_t +bsp_construct_array_t_allocator(bsp_array_t* array, size_t size, + bsp_type_t type, bsp_allocator_t allocator) { size_t byte_size = size * bsp_type_size(type); - bsp_array_t array; - array.allocator = bsp_default_allocator; - array.data = array.allocator.malloc(byte_size); - assert(array.data != NULL); - array.size = size; - array.type = type; + array->allocator = allocator; + array->data = array->allocator.malloc(byte_size); + + if (array->data == NULL) { + return BSP_ERROR_MEMORY; + } - return array; + array->size = size; + array->type = type; + + return BSP_SUCCESS; } -static inline bsp_array_t bsp_copy_construct_array_t(bsp_array_t other) { - bsp_array_t array = bsp_construct_array_t(other.size, other.type); - memcpy(array.data, other.data, other.size * bsp_type_size(other.type)); +static inline bsp_error_t bsp_construct_array_t(bsp_array_t* array, size_t size, + bsp_type_t type) { + return bsp_construct_array_t_allocator(array, size, type, + bsp_default_allocator); +} - return array; +static inline bsp_error_t bsp_copy_construct_array_t(bsp_array_t* array, + bsp_array_t other) { + bsp_error_t error = bsp_construct_array_t(array, other.size, other.type); + if (error != BSP_SUCCESS) { + return error; + } + + memcpy(array->data, other.data, other.size * bsp_type_size(other.type)); + + return BSP_SUCCESS; } -static inline bsp_array_t bsp_complex_array_to_fp(bsp_array_t other) { - assert(other.type == BSP_COMPLEX_FLOAT32 || - other.type == BSP_COMPLEX_FLOAT64); +static inline bsp_error_t bsp_complex_array_to_fp(bsp_array_t* array) { + if (array->type != BSP_COMPLEX_FLOAT32 && + array->type != BSP_COMPLEX_FLOAT64) { + return BSP_ERROR_TYPE; + } - bsp_array_t array; - array.data = other.data; - array.size = other.size * 2; - array.allocator = other.allocator; + array->size = array->size * 2; - if (other.type == BSP_COMPLEX_FLOAT32) { - array.type = BSP_FLOAT32; + if (array->type == BSP_COMPLEX_FLOAT32) { + array->type = BSP_FLOAT32; } else { - array.type = BSP_FLOAT64; + array->type = BSP_FLOAT64; } - return array; + return BSP_SUCCESS; } -static inline bsp_array_t bsp_fp_array_to_complex(bsp_array_t other) { - assert(other.type == BSP_FLOAT32 || other.type == BSP_FLOAT64); - - bsp_array_t array; - array.data = other.data; - array.size = other.size / 2; - array.allocator = other.allocator; +static inline bsp_error_t bsp_fp_array_to_complex(bsp_array_t* array) { + if (array->type != BSP_FLOAT32 && array->type != BSP_FLOAT64) { + return BSP_ERROR_TYPE; + } - if (other.type == BSP_FLOAT32) { - array.type = BSP_COMPLEX_FLOAT32; + if (array->type == BSP_FLOAT32) { + array->type = BSP_COMPLEX_FLOAT32; } else { - array.type = BSP_COMPLEX_FLOAT64; + array->type = BSP_COMPLEX_FLOAT64; } - return array; + array->size = array->size / 2; + + return BSP_SUCCESS; } -static inline void bsp_destroy_array_t(bsp_array_t array) { - array.allocator.free(array.data); +static inline void bsp_destroy_array_t(bsp_array_t* array) { + array->allocator.free(array->data); } static inline bool bsp_array_equal(bsp_array_t x, bsp_array_t y) { @@ -145,7 +165,7 @@ static inline bool bsp_array_equal(bsp_array_t x, bsp_array_t y) { float* ptr = (float*) array.data; \ code_block; \ } else if (array.type == BSP_FLOAT64) { \ - double* ptr = (float*) array.data; \ + double* ptr = (double*) array.data; \ code_block; \ } else if (array.type == BSP_BINT8) { \ int8_t* ptr = (int8_t*) array.data; \ diff --git a/include/binsparse/binsparse.h b/include/binsparse/binsparse.h index bb6eecd..189aa7c 100644 --- a/include/binsparse/binsparse.h +++ b/include/binsparse/binsparse.h @@ -11,8 +11,8 @@ #include #include #include +#include #include -#include #include #include #include diff --git a/include/binsparse/binsparse_all.h b/include/binsparse/binsparse_all.h new file mode 100644 index 0000000..34f3463 --- /dev/null +++ b/include/binsparse/binsparse_all.h @@ -0,0 +1,11 @@ +/* + * SPDX-FileCopyrightText: 2024 Binsparse Developers + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#pragma once + +#include +#include +#include diff --git a/include/binsparse/binsparse_cJSON.h b/include/binsparse/binsparse_cJSON.h new file mode 100644 index 0000000..c74b9c4 --- /dev/null +++ b/include/binsparse/binsparse_cJSON.h @@ -0,0 +1,34 @@ +/* + * SPDX-FileCopyrightText: 2024 Binsparse Developers + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#pragma once + +#include +#include + +#include + +#ifndef BSP_BINSPARSE_CJSON_H +#define BSP_BINSPARSE_CJSON_H +#endif + +#ifdef __cplusplus +extern "C" { +#endif + +bsp_error_t bsp_write_matrix_cjson(const char* fname, bsp_matrix_t matrix, + const char* group, cJSON* user_json, + int compression_level); + +#ifdef BSP_BINSPARSE_HDF5_H +bsp_error_t bsp_write_matrix_to_group_cjson(hid_t f, bsp_matrix_t matrix, + cJSON* user_json, + int compression_level); +#endif + +#ifdef __cplusplus +} +#endif diff --git a/include/binsparse/binsparse_hdf5.h b/include/binsparse/binsparse_hdf5.h new file mode 100644 index 0000000..2aef4a5 --- /dev/null +++ b/include/binsparse/binsparse_hdf5.h @@ -0,0 +1,45 @@ +/* + * SPDX-FileCopyrightText: 2024 Binsparse Developers + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#pragma once + +#include +#include + +#include + +#include + +#ifndef BSP_BINSPARSE_HDF5_H +#define BSP_BINSPARSE_HDF5_H +#endif + +#ifdef __cplusplus +extern "C" { +#endif + +#if __STDC_VERSION__ >= 201112L +bsp_error_t bsp_read_matrix_from_group_parallel(bsp_matrix_t* matrix, hid_t f, + int num_threads); +#endif + +bsp_error_t bsp_read_matrix_from_group(bsp_matrix_t* matrix, hid_t f); +bsp_error_t bsp_read_matrix_from_group_allocator(bsp_matrix_t* matrix, hid_t f, + bsp_allocator_t allocator); + +bsp_error_t bsp_write_matrix_to_group(hid_t f, bsp_matrix_t matrix, + const char* user_json, + int compression_level); + +#ifdef BSP_BINSPARSE_CJSON_H +bsp_error_t bsp_write_matrix_to_group_cjson(hid_t f, bsp_matrix_t matrix, + cJSON* user_json, + int compression_level); +#endif + +#ifdef __cplusplus +} +#endif diff --git a/include/binsparse/convert_matrix.h b/include/binsparse/convert_matrix.h index 2f6c332..3d2c12d 100644 --- a/include/binsparse/convert_matrix.h +++ b/include/binsparse/convert_matrix.h @@ -8,9 +8,26 @@ #include #include +#include +#include +#include + +static inline bsp_error_t +bsp_copy_construct_array_t_allocator(bsp_array_t* array, bsp_array_t other, + bsp_allocator_t allocator) { + bsp_error_t error = + bsp_construct_array_t_allocator(array, other.size, other.type, allocator); + if (error != BSP_SUCCESS) { + return error; + } -static inline bsp_matrix_t bsp_convert_matrix(bsp_matrix_t matrix, - bsp_matrix_format_t format) { + memcpy(array->data, other.data, other.size * bsp_type_size(other.type)); + return BSP_SUCCESS; +} + +static inline bsp_matrix_t +bsp_convert_matrix_allocator(bsp_matrix_t matrix, bsp_matrix_format_t format, + bsp_allocator_t allocator) { // Throw an error if matrix already in desired format. if (matrix.format == format) { assert(false); @@ -20,7 +37,8 @@ static inline bsp_matrix_t bsp_convert_matrix(bsp_matrix_t matrix, // *Convert to COO* from another format. if (matrix.format == BSP_CSR) { // Convert CSR -> COOR - bsp_matrix_t result = bsp_construct_default_matrix_t(); + bsp_matrix_t result; + bsp_construct_default_matrix_t_allocator(&result, allocator); result.format = BSP_COOR; @@ -37,17 +55,37 @@ static inline bsp_matrix_t bsp_convert_matrix(bsp_matrix_t matrix, bsp_type_t index_type = bsp_pick_integer_type(max_dim); - result.values = bsp_copy_construct_array_t(matrix.values); + bsp_error_t error = bsp_copy_construct_array_t_allocator( + &result.values, matrix.values, allocator); + if (error != BSP_SUCCESS) { + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } // There is a corner case with tall and skinny matrices where we need a // higher width for rowind. In order to keep rowind/colind the same type, // we might upcast. if (index_type == matrix.indices_1.type) { - result.indices_1 = bsp_copy_construct_array_t(matrix.indices_1); + error = bsp_copy_construct_array_t_allocator( + &result.indices_1, matrix.indices_1, allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&result.values); + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } } else { - result.indices_1 = - bsp_construct_array_t(matrix.indices_1.size, index_type); + error = bsp_construct_array_t_allocator( + &result.indices_1, matrix.indices_1.size, index_type, allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&result.values); + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } + for (size_t i = 0; i < matrix.indices_1.size; i++) { size_t index; bsp_array_read(matrix.indices_1, i, index); @@ -55,7 +93,15 @@ static inline bsp_matrix_t bsp_convert_matrix(bsp_matrix_t matrix, } } - result.indices_0 = bsp_construct_array_t(matrix.nnz, index_type); + error = bsp_construct_array_t_allocator(&result.indices_0, matrix.nnz, + index_type, allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&result.values); + bsp_destroy_array_t(&result.indices_1); + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } for (size_t i = 0; i < matrix.nrows; i++) { size_t row_begin, row_end; @@ -66,6 +112,168 @@ static inline bsp_matrix_t bsp_convert_matrix(bsp_matrix_t matrix, } } return result; + } else if (matrix.format == BSP_CSC) { + // Convert CSC -> COOR + bsp_matrix_t result; + bsp_construct_default_matrix_t_allocator(&result, allocator); + + result.format = BSP_COOR; + + // Inherit NNZ, nrows, ncols, ISO-ness, and structure directly from + // original matrix. + result.nnz = matrix.nnz; + result.nrows = matrix.nrows; + result.ncols = matrix.ncols; + result.is_iso = matrix.is_iso; + result.structure = matrix.structure; + + size_t max_dim = + (matrix.nrows > matrix.ncols) ? matrix.nrows : matrix.ncols; + + bsp_type_t index_type = bsp_pick_integer_type(max_dim); + + bsp_error_t error = bsp_copy_construct_array_t_allocator( + &result.values, matrix.values, allocator); + if (error != BSP_SUCCESS) { + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } + + // Copy row indices from CSC to become row indices in COOR. + if (index_type == matrix.indices_1.type) { + error = bsp_copy_construct_array_t_allocator( + &result.indices_0, matrix.indices_1, allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&result.values); + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } + } else { + error = bsp_construct_array_t_allocator( + &result.indices_0, matrix.indices_1.size, index_type, allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&result.values); + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } + + for (size_t i = 0; i < matrix.indices_1.size; i++) { + size_t index; + bsp_array_read(matrix.indices_1, i, index); + bsp_array_write(result.indices_0, i, index); + } + } + + // Generate column indices by expanding column pointers. + error = bsp_construct_array_t_allocator(&result.indices_1, matrix.nnz, + index_type, allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&result.values); + bsp_destroy_array_t(&result.indices_0); + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } + + for (size_t j = 0; j < matrix.ncols; j++) { + size_t col_begin, col_end; + bsp_array_read(matrix.pointers_to_1, j, col_begin); + bsp_array_read(matrix.pointers_to_1, j + 1, col_end); + for (size_t i_ptr = col_begin; i_ptr < col_end; i_ptr++) { + bsp_array_write(result.indices_1, i_ptr, j); + } + } + + // Sort the result by rows to produce valid COOR. + size_t* indices = (size_t*) malloc(sizeof(size_t) * matrix.nnz); + if (indices == NULL) { + bsp_destroy_array_t(&result.values); + bsp_destroy_array_t(&result.indices_0); + bsp_destroy_array_t(&result.indices_1); + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } + + for (size_t i = 0; i < matrix.nnz; i++) { + indices[i] = i; + } + + bsp_coo_indices_.rowind = result.indices_0; + bsp_coo_indices_.colind = result.indices_1; + + qsort(indices, matrix.nnz, sizeof(size_t), + bsp_coo_comparison_row_sort_operator_impl_); + + bsp_array_t rowind; + bsp_array_t colind; + + error = bsp_copy_construct_array_t_allocator(&rowind, result.indices_0, + allocator); + if (error != BSP_SUCCESS) { + free(indices); + bsp_destroy_array_t(&result.values); + bsp_destroy_array_t(&result.indices_0); + bsp_destroy_array_t(&result.indices_1); + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } + + error = bsp_copy_construct_array_t_allocator(&colind, result.indices_1, + allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&rowind); + free(indices); + bsp_destroy_array_t(&result.values); + bsp_destroy_array_t(&result.indices_0); + bsp_destroy_array_t(&result.indices_1); + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } + + bsp_array_t values; + + if (!result.is_iso) { + error = bsp_copy_construct_array_t_allocator(&values, result.values, + allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&rowind); + bsp_destroy_array_t(&colind); + free(indices); + bsp_destroy_array_t(&result.values); + bsp_destroy_array_t(&result.indices_0); + bsp_destroy_array_t(&result.indices_1); + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } + } + + for (size_t i = 0; i < matrix.nnz; i++) { + bsp_array_awrite(rowind, i, result.indices_0, indices[i]); + bsp_array_awrite(colind, i, result.indices_1, indices[i]); + if (!result.is_iso) { + bsp_array_awrite(values, i, result.values, indices[i]); + } + } + + bsp_destroy_array_t(&result.indices_0); + bsp_destroy_array_t(&result.indices_1); + result.indices_0 = rowind; + result.indices_1 = colind; + + if (!result.is_iso) { + bsp_destroy_array_t(&result.values); + result.values = values; + } + + free(indices); + return result; } else { assert(false); } @@ -75,15 +283,18 @@ static inline bsp_matrix_t bsp_convert_matrix(bsp_matrix_t matrix, // Currently only support COOR -> X. // If matrix is not COOR, convert to COOR. if (matrix.format != BSP_COOR) { - bsp_matrix_t intermediate = bsp_convert_matrix(matrix, BSP_COOR); - bsp_matrix_t result = bsp_convert_matrix(intermediate, format); - bsp_destroy_matrix_t(intermediate); + bsp_matrix_t intermediate = + bsp_convert_matrix_allocator(matrix, BSP_COOR, allocator); + bsp_matrix_t result = + bsp_convert_matrix_allocator(intermediate, format, allocator); + bsp_destroy_matrix_t(&intermediate); return result; } else { if (format == BSP_CSR) { // Convert COOR -> CSR - bsp_matrix_t result = bsp_construct_default_matrix_t(); + bsp_matrix_t result; + bsp_construct_default_matrix_t_allocator(&result, allocator); result.format = BSP_CSR; @@ -109,12 +320,32 @@ static inline bsp_matrix_t bsp_convert_matrix(bsp_matrix_t matrix, // indices can be copied exactly. Values' type will not change, but // column indices might, thus the extra branch. - result.values = bsp_copy_construct_array_t(matrix.values); + bsp_error_t error = bsp_copy_construct_array_t_allocator( + &result.values, matrix.values, allocator); + if (error != BSP_SUCCESS) { + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } if (index_type == matrix.indices_1.type) { - result.indices_1 = bsp_copy_construct_array_t(matrix.indices_1); + error = bsp_copy_construct_array_t_allocator( + &result.indices_1, matrix.indices_1, allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&result.values); + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } } else { - result.indices_1 = bsp_construct_array_t(matrix.nnz, index_type); + error = bsp_construct_array_t_allocator(&result.indices_1, matrix.nnz, + index_type, allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&result.values); + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } for (size_t i = 0; i < matrix.nnz; i++) { size_t index; @@ -123,8 +354,15 @@ static inline bsp_matrix_t bsp_convert_matrix(bsp_matrix_t matrix, } } - result.pointers_to_1 = - bsp_construct_array_t(matrix.nrows + 1, index_type); + error = bsp_construct_array_t_allocator( + &result.pointers_to_1, matrix.nrows + 1, index_type, allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&result.values); + bsp_destroy_array_t(&result.indices_1); + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } bsp_array_t rowptr = result.pointers_to_1; @@ -148,6 +386,111 @@ static inline bsp_matrix_t bsp_convert_matrix(bsp_matrix_t matrix, bsp_array_write(rowptr, r + 1, matrix.nnz); } + return result; + } else if (format == BSP_CSC) { + // Convert COOR -> CSC + + // First, sort by columns to prepare for CSC format. + size_t* indices = (size_t*) malloc(sizeof(size_t) * matrix.nnz); + if (indices == NULL) { + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } + + for (size_t i = 0; i < matrix.nnz; i++) { + indices[i] = i; + } + + bsp_coo_indices_.rowind = matrix.indices_0; + bsp_coo_indices_.colind = matrix.indices_1; + + qsort(indices, matrix.nnz, sizeof(size_t), + bsp_coo_comparison_col_sort_operator_impl_); + + bsp_matrix_t result; + bsp_construct_default_matrix_t_allocator(&result, allocator); + + result.format = BSP_CSC; + + result.nrows = matrix.nrows; + result.ncols = matrix.ncols; + result.nnz = matrix.nnz; + result.is_iso = matrix.is_iso; + result.structure = matrix.structure; + + size_t max_dim = + (matrix.nrows > matrix.ncols) ? matrix.nrows : matrix.ncols; + + size_t max_value = + (max_dim > matrix.values.size) ? max_dim : matrix.values.size; + + bsp_type_t index_type = bsp_pick_integer_type(max_value); + + // Reorder values according to column-major sort. + bsp_error_t error = bsp_construct_array_t_allocator( + &result.values, matrix.values.size, matrix.values.type, allocator); + if (error != BSP_SUCCESS) { + free(indices); + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } + + for (size_t i = 0; i < matrix.values.size; i++) { + bsp_array_awrite(result.values, i, matrix.values, + matrix.is_iso ? 0 : indices[i]); + } + + // Reorder row indices according to column-major sort. + error = bsp_construct_array_t_allocator(&result.indices_1, matrix.nnz, + index_type, allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&result.values); + free(indices); + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } + + for (size_t i = 0; i < matrix.nnz; i++) { + bsp_array_awrite(result.indices_1, i, matrix.indices_0, indices[i]); + } + + // Build column pointers. + error = bsp_construct_array_t_allocator( + &result.pointers_to_1, matrix.ncols + 1, index_type, allocator); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&result.values); + bsp_destroy_array_t(&result.indices_1); + free(indices); + bsp_matrix_t empty_result; + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); + return empty_result; + } + + bsp_array_t colptr = result.pointers_to_1; + + bsp_array_write(colptr, 0, 0); + + size_t c = 0; + for (size_t i = 0; i < matrix.nnz; i++) { + size_t j; + bsp_array_read(matrix.indices_1, indices[i], j); + + while (c < j) { + assert(c + 1 <= matrix.ncols); + + bsp_array_write(colptr, c + 1, i); + c++; + } + } + + for (; c < matrix.ncols; c++) { + bsp_array_write(colptr, c + 1, matrix.nnz); + } + + free(indices); return result; } else { assert(false); @@ -155,3 +498,8 @@ static inline bsp_matrix_t bsp_convert_matrix(bsp_matrix_t matrix, } } } + +static inline bsp_matrix_t bsp_convert_matrix(bsp_matrix_t matrix, + bsp_matrix_format_t format) { + return bsp_convert_matrix_allocator(matrix, format, bsp_default_allocator); +} diff --git a/include/binsparse/detail/hdf5_types.h b/include/binsparse/detail/hdf5_types.h new file mode 100644 index 0000000..02f0315 --- /dev/null +++ b/include/binsparse/detail/hdf5_types.h @@ -0,0 +1,208 @@ +/* + * SPDX-FileCopyrightText: 2024 Binsparse Developers + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#pragma once + +#include +#include + +static inline hid_t bsp_get_hdf5_standard_type(bsp_type_t type) { + if (type == BSP_UINT8) { + return H5T_STD_U8LE; + } else if (type == BSP_UINT16) { + return H5T_STD_U16LE; + } else if (type == BSP_UINT32) { + return H5T_STD_U32LE; + } else if (type == BSP_UINT64) { + return H5T_STD_U64LE; + } else if (type == BSP_INT8) { + return H5T_STD_I8LE; + } else if (type == BSP_INT16) { + return H5T_STD_I16LE; + } else if (type == BSP_INT32) { + return H5T_STD_I32LE; + } else if (type == BSP_INT64) { + return H5T_STD_I64LE; + } else if (type == BSP_FLOAT32) { + return H5T_IEEE_F32LE; + } else if (type == BSP_FLOAT64) { + return H5T_IEEE_F64LE; + } else if (type == BSP_BINT8) { + return H5T_STD_I8LE; + } else { + return H5I_INVALID_HID; + } +} + +static inline bsp_type_t bsp_get_bsp_type(hid_t type) { + H5T_class_t cl = H5Tget_class(type); + H5T_order_t order = H5Tget_order(type); + H5T_sign_t sign = H5Tget_sign(type); + size_t size = H5Tget_size(type); + + if (cl == H5T_INTEGER) { + if (sign == H5T_SGN_NONE) { + if (size == 1) { + return BSP_UINT8; + } else if (size == 2) { + return BSP_UINT16; + } else if (size == 4) { + return BSP_UINT32; + } else if (size == 8) { + return BSP_UINT64; + } else { + return BSP_INVALID_TYPE; + } + } else /* if (sign == H5T_SGN_2) */ { + if (size == 1) { + return BSP_INT8; + } else if (size == 2) { + return BSP_INT16; + } else if (size == 4) { + return BSP_INT32; + } else if (size == 8) { + return BSP_INT64; + } else { + return BSP_INVALID_TYPE; + } + } + } else if (cl == H5T_FLOAT) { + if (size == 4) { + return BSP_FLOAT32; + } else if (size == 8) { + return BSP_FLOAT64; + } else { + return BSP_INVALID_TYPE; + } + } else { + return BSP_INVALID_TYPE; + } +} + +// NOTE: This code is a bit silly, but it seems to be the only +// way to generically determine the HDF5 native types for +// stdint's fixed width integer types. +static inline hid_t bsp_get_hdf5_native_type(bsp_type_t type) { + if (type == BSP_INT8 || type == BSP_BINT8) { + if (sizeof(int8_t) == sizeof(char)) { + return H5T_NATIVE_CHAR; + } else if (sizeof(int8_t) == sizeof(short)) { + return H5T_NATIVE_SHORT; + } else if (sizeof(int8_t) == sizeof(int)) { + return H5T_NATIVE_INT; + } else if (sizeof(int8_t) == sizeof(long)) { + return H5T_NATIVE_LONG; + } else if (sizeof(int8_t) == sizeof(long long)) { + return H5T_NATIVE_LLONG; + } else { + assert(false); + } + } else if (type == BSP_INT16) { + if (sizeof(int16_t) == sizeof(char)) { + return H5T_NATIVE_CHAR; + } else if (sizeof(int16_t) == sizeof(short)) { + return H5T_NATIVE_SHORT; + } else if (sizeof(int16_t) == sizeof(int)) { + return H5T_NATIVE_INT; + } else if (sizeof(int32_t) == sizeof(long)) { + return H5T_NATIVE_LONG; + } else if (sizeof(int64_t) == sizeof(long long)) { + return H5T_NATIVE_LLONG; + } else { + assert(false); + } + } else if (type == BSP_INT32) { + if (sizeof(int32_t) == sizeof(char)) { + return H5T_NATIVE_CHAR; + } else if (sizeof(int32_t) == sizeof(short)) { + return H5T_NATIVE_SHORT; + } else if (sizeof(int32_t) == sizeof(int)) { + return H5T_NATIVE_INT; + } else if (sizeof(int32_t) == sizeof(long)) { + return H5T_NATIVE_LONG; + } else if (sizeof(int32_t) == sizeof(long long)) { + return H5T_NATIVE_LLONG; + } else { + assert(false); + } + } else if (type == BSP_INT64) { + if (sizeof(int64_t) == sizeof(char)) { + return H5T_NATIVE_CHAR; + } else if (sizeof(int64_t) == sizeof(short)) { + return H5T_NATIVE_SHORT; + } else if (sizeof(int64_t) == sizeof(int)) { + return H5T_NATIVE_INT; + } else if (sizeof(int64_t) == sizeof(long)) { + return H5T_NATIVE_LONG; + } else if (sizeof(int64_t) == sizeof(long long)) { + return H5T_NATIVE_LLONG; + } else { + assert(false); + } + } else if (type == BSP_UINT8) { + if (sizeof(uint8_t) == sizeof(unsigned char)) { + return H5T_NATIVE_UCHAR; + } else if (sizeof(uint8_t) == sizeof(unsigned short)) { + return H5T_NATIVE_USHORT; + } else if (sizeof(uint8_t) == sizeof(unsigned int)) { + return H5T_NATIVE_UINT; + } else if (sizeof(uint8_t) == sizeof(unsigned long)) { + return H5T_NATIVE_ULONG; + } else if (sizeof(uint8_t) == sizeof(unsigned long long)) { + return H5T_NATIVE_ULLONG; + } else { + assert(false); + } + } else if (type == BSP_UINT16) { + if (sizeof(uint16_t) == sizeof(unsigned char)) { + return H5T_NATIVE_UCHAR; + } else if (sizeof(uint16_t) == sizeof(unsigned short)) { + return H5T_NATIVE_USHORT; + } else if (sizeof(uint16_t) == sizeof(unsigned int)) { + return H5T_NATIVE_UINT; + } else if (sizeof(uint16_t) == sizeof(unsigned long)) { + return H5T_NATIVE_ULONG; + } else if (sizeof(uint16_t) == sizeof(unsigned long long)) { + return H5T_NATIVE_ULLONG; + } else { + assert(false); + } + } else if (type == BSP_UINT32) { + if (sizeof(uint32_t) == sizeof(unsigned char)) { + return H5T_NATIVE_UCHAR; + } else if (sizeof(uint32_t) == sizeof(unsigned short)) { + return H5T_NATIVE_USHORT; + } else if (sizeof(uint32_t) == sizeof(unsigned int)) { + return H5T_NATIVE_UINT; + } else if (sizeof(uint32_t) == sizeof(unsigned long)) { + return H5T_NATIVE_ULONG; + } else if (sizeof(uint32_t) == sizeof(unsigned long long)) { + return H5T_NATIVE_ULLONG; + } else { + assert(false); + } + } else if (type == BSP_UINT64) { + if (sizeof(uint64_t) == sizeof(unsigned char)) { + return H5T_NATIVE_UCHAR; + } else if (sizeof(uint64_t) == sizeof(unsigned short)) { + return H5T_NATIVE_USHORT; + } else if (sizeof(uint64_t) == sizeof(unsigned int)) { + return H5T_NATIVE_UINT; + } else if (sizeof(uint64_t) == sizeof(unsigned long)) { + return H5T_NATIVE_ULONG; + } else if (sizeof(uint64_t) == sizeof(unsigned long long)) { + return H5T_NATIVE_ULLONG; + } else { + assert(false); + } + } else if (type == BSP_FLOAT32) { + return H5T_NATIVE_FLOAT; + } else if (type == BSP_FLOAT64) { + return H5T_NATIVE_DOUBLE; + } else { + return H5I_INVALID_HID; + } +} diff --git a/include/binsparse/detail/parse_dataset.h b/include/binsparse/detail/parse_dataset.h index 4facff7..489be88 100644 --- a/include/binsparse/detail/parse_dataset.h +++ b/include/binsparse/detail/parse_dataset.h @@ -13,6 +13,11 @@ typedef struct { char* dataset; } bsp_fdataset_info_t; +static inline void bsp_destroy_fdataset_info_t(bsp_fdataset_info_t* info) { + free(info->fname); + free(info->dataset); +} + static inline bsp_fdataset_info_t bsp_parse_fdataset_string(char* str) { size_t len = strlen(str); diff --git a/include/binsparse/error.h b/include/binsparse/error.h new file mode 100644 index 0000000..73b3f92 --- /dev/null +++ b/include/binsparse/error.h @@ -0,0 +1,83 @@ +/* + * SPDX-FileCopyrightText: 2024 Binsparse Developers + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +#pragma once + +#include +#include + +typedef enum bsp_error_t { + BSP_SUCCESS = 0, + + // Memory-related failures (malloc, array allocation, etc.) + BSP_ERROR_MEMORY = 1, + + // File I/O failures (file not found, read/write errors, HDF5 operations) + BSP_ERROR_IO = 2, + + // JSON parsing and metadata errors + BSP_ERROR_FORMAT = 3, + + // Invalid input parameters, dimensions, types, etc. + BSP_ERROR_INVALID_INPUT = 4, + + // Unsupported operations or not-yet-implemented features + BSP_ERROR_UNSUPPORTED = 5, + + // Type-related errors (invalid types, type mismatches) + BSP_ERROR_TYPE = 6, + + // Generic internal errors (should-never-happen cases) + BSP_ERROR_INTERNAL = 7 + +} bsp_error_t; + +/** + * Get a human-readable error message for a bsp_error_t code. + * + * @param error The error code + * @return A constant string describing the error + */ +static inline const char* bsp_get_error_string(bsp_error_t error) { + switch (error) { + case BSP_SUCCESS: + return "Success"; + case BSP_ERROR_MEMORY: + return "Memory allocation or management error"; + case BSP_ERROR_IO: + return "File I/O or HDF5 operation error"; + case BSP_ERROR_FORMAT: + return "JSON parsing or metadata format error"; + case BSP_ERROR_INVALID_INPUT: + return "Invalid input parameters or data"; + case BSP_ERROR_UNSUPPORTED: + return "Unsupported operation or feature"; + case BSP_ERROR_TYPE: + return "Data type error or type mismatch"; + case BSP_ERROR_INTERNAL: + return "Internal library error"; + default: + return "Unknown error"; + } +} + +/** + * BSP_CHECK macro - checks if a bsp_error_t is BSP_SUCCESS. + * If not, prints file, line number, and error message, then aborts. + * + * Usage: BSP_CHECK(bsp_read_matrix(&matrix, "file.hdf5", NULL)); + * + * @param call The function call that returns a bsp_error_t + */ +#define BSP_CHECK(call) \ + do { \ + bsp_error_t _bsp_error = (call); \ + if (_bsp_error != BSP_SUCCESS) { \ + fprintf(stderr, "BSP Error at %s:%d: %s\n", __FILE__, __LINE__, \ + bsp_get_error_string(_bsp_error)); \ + abort(); \ + } \ + } while (0) diff --git a/include/binsparse/generate.h b/include/binsparse/generate.h index 93205a1..b625669 100644 --- a/include/binsparse/generate.h +++ b/include/binsparse/generate.h @@ -70,13 +70,31 @@ static inline void bsp_array_fill_random(bsp_array_t array, size_t bound) { static inline bsp_matrix_t bsp_generate_coo(size_t m, size_t n, size_t nnz, bsp_type_t value_type, bsp_type_t index_type) { - bsp_matrix_t matrix = bsp_construct_default_matrix_t(); + bsp_matrix_t matrix; + bsp_construct_default_matrix_t(&matrix); matrix.nrows = m; matrix.ncols = n; matrix.nnz = nnz; - matrix.values = bsp_construct_array_t(nnz, value_type); - matrix.indices_0 = bsp_construct_array_t(nnz, index_type); - matrix.indices_1 = bsp_construct_array_t(nnz, index_type); + + bsp_error_t error = bsp_construct_array_t(&matrix.values, nnz, value_type); + if (error != BSP_SUCCESS) { + // Return empty matrix on error + return matrix; + } + + error = bsp_construct_array_t(&matrix.indices_0, nnz, index_type); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&matrix.values); + return matrix; + } + + error = bsp_construct_array_t(&matrix.indices_1, nnz, index_type); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&matrix.values); + bsp_destroy_array_t(&matrix.indices_0); + return matrix; + } + matrix.format = BSP_COO; bsp_array_fill_random(matrix.values, 100); diff --git a/include/binsparse/hdf5_wrapper.h b/include/binsparse/hdf5_wrapper.h index 5aded28..f8b44df 100644 --- a/include/binsparse/hdf5_wrapper.h +++ b/include/binsparse/hdf5_wrapper.h @@ -13,6 +13,7 @@ #include +#include #include #if __STDC_VERSION__ >= 201112L @@ -21,10 +22,14 @@ // Write an array to a dataset / file // Returns 0 on success, nonzero on error. -static inline int bsp_write_array(hid_t f, const char* label, bsp_array_t array, - int compression_level) { +static inline bsp_error_t bsp_write_array(hid_t f, const char* label, + bsp_array_t array, + int compression_level) { if (array.type == BSP_COMPLEX_FLOAT32 || array.type == BSP_COMPLEX_FLOAT64) { - array = bsp_complex_array_to_fp(array); + bsp_error_t error = bsp_complex_array_to_fp(&array); + if (error != BSP_SUCCESS) { + return BSP_ERROR_TYPE; + } } hsize_t hsize[1]; @@ -56,7 +61,10 @@ static inline int bsp_write_array(hid_t f, const char* label, bsp_array_t array, H5Dcreate2(f, label, hdf5_standard_type, fspace, lcpl, dcpl, H5P_DEFAULT); if (dset == H5I_INVALID_HID) { - return -1; + H5Sclose(fspace); + H5Pclose(lcpl); + H5Pclose(dcpl); + return BSP_ERROR_IO; } hid_t hdf5_native_type = bsp_get_hdf5_native_type(array.type); @@ -64,30 +72,41 @@ static inline int bsp_write_array(hid_t f, const char* label, bsp_array_t array, hid_t r = H5Dwrite(dset, hdf5_native_type, H5S_ALL, fspace, H5P_DEFAULT, array.data); - if (r == H5I_INVALID_HID) { - return -2; + if (r < 0) { + H5Dclose(dset); + H5Sclose(fspace); + H5Pclose(lcpl); + H5Pclose(dcpl); + return BSP_ERROR_IO; } H5Sclose(fspace); + H5Dclose(dset); H5Pclose(lcpl); H5Pclose(dcpl); - return 0; + return BSP_SUCCESS; } #if __STDC_VERSION__ >= 201112L -static inline bsp_array_t bsp_read_array_parallel(hid_t f, const char* label, +static inline bsp_error_t bsp_read_array_parallel(bsp_array_t* array, hid_t f, + const char* label, int num_threads) { hid_t dset = H5Dopen2(f, label, H5P_DEFAULT); if (dset == H5I_INVALID_HID) { - return bsp_construct_default_array_t(); + H5Dclose(dset); + bsp_construct_default_array_t(array); + return BSP_ERROR_IO; } hid_t fspace = H5Dget_space(dset); if (fspace == H5I_INVALID_HID) { - return bsp_construct_default_array_t(); + H5Sclose(fspace); + H5Dclose(dset); + bsp_construct_default_array_t(array); + return BSP_ERROR_IO; } hsize_t dims[3]; @@ -95,7 +114,10 @@ static inline bsp_array_t bsp_read_array_parallel(hid_t f, const char* label, int r = H5Sget_simple_extent_dims(fspace, dims, NULL); if (r < 0) { - return bsp_construct_default_array_t(); + H5Sclose(fspace); + H5Dclose(dset); + bsp_construct_default_array_t(array); + return BSP_ERROR_IO; } hid_t hdf5_type = H5Dget_type(dset); @@ -104,10 +126,9 @@ static inline bsp_array_t bsp_read_array_parallel(hid_t f, const char* label, // Array will be written into a POSIX shared memory. bsp_shm_t array_shm = bsp_shm_new(dims[0] * bsp_type_size(type)); - bsp_array_t array; - array.type = type; - array.size = dims[0]; - array.allocator = bsp_shm_allocator; + array->type = type; + array->size = dims[0]; + array->allocator = bsp_shm_allocator; bsp_shm_t active_children_shm = bsp_shm_new(sizeof(_Atomic int)); @@ -131,17 +152,17 @@ static inline bsp_array_t bsp_read_array_parallel(hid_t f, const char* label, } } - array.data = bsp_shm_attach(array_shm); + array->data = bsp_shm_attach(array_shm); if (thread_num == 0) { bsp_shm_delete(array_shm); } - hsize_t chunk_size = (array.size + num_threads - 1) / num_threads; - hsize_t start = (chunk_size * thread_num < array.size) + hsize_t chunk_size = (array->size + num_threads - 1) / num_threads; + hsize_t start = (chunk_size * thread_num < array->size) ? chunk_size * thread_num - : array.size; + : array->size; hsize_t count = - (start + chunk_size <= array.size) ? chunk_size : array.size - start; + (start + chunk_size <= array->size) ? chunk_size : array->size - start; if (count > 0) { H5Sselect_hyperslab(fspace, H5S_SELECT_SET, &start, NULL, &count, NULL); @@ -149,17 +170,17 @@ static inline bsp_array_t bsp_read_array_parallel(hid_t f, const char* label, hid_t memspace_id = H5Screate_simple(1, &count, NULL); H5Dread(dset, bsp_get_hdf5_native_type(type), memspace_id, fspace, - H5P_DEFAULT, ((char*) array.data) + start * bsp_type_size(type)); + H5P_DEFAULT, ((char*) array->data) + start * bsp_type_size(type)); H5Sclose(memspace_id); } - H5Dclose(dset); H5Sclose(fspace); + H5Dclose(dset); if (thread_num > 0) { atomic_fetch_add_explicit(active_children, -1, memory_order_relaxed); bsp_shm_detach(active_children); - bsp_shm_detach(array.data); + bsp_shm_detach(array->data); exit(0); } @@ -169,21 +190,27 @@ static inline bsp_array_t bsp_read_array_parallel(hid_t f, const char* label, } bsp_shm_detach(active_children); - return array; + return BSP_SUCCESS; } #endif -static inline bsp_array_t bsp_read_array(hid_t f, const char* label) { +static inline bsp_error_t bsp_read_array_allocator(bsp_array_t* array, hid_t f, + const char* label, + bsp_allocator_t allocator) { hid_t dset = H5Dopen2(f, label, H5P_DEFAULT); if (dset == H5I_INVALID_HID) { - return bsp_construct_default_array_t(); + bsp_construct_default_array_t_allocator(array, allocator); + return BSP_ERROR_IO; } hid_t fspace = H5Dget_space(dset); if (fspace == H5I_INVALID_HID) { - return bsp_construct_default_array_t(); + H5Sclose(fspace); + H5Dclose(dset); + bsp_construct_default_array_t_allocator(array, allocator); + return BSP_ERROR_IO; } hsize_t dims[3]; @@ -191,29 +218,48 @@ static inline bsp_array_t bsp_read_array(hid_t f, const char* label) { int r = H5Sget_simple_extent_dims(fspace, dims, NULL); if (r < 0) { - return bsp_construct_default_array_t(); + H5Dclose(dset); + H5Sclose(fspace); + bsp_construct_default_array_t_allocator(array, allocator); + return BSP_ERROR_IO; } hid_t hdf5_type = H5Dget_type(dset); bsp_type_t type = bsp_get_bsp_type(hdf5_type); - bsp_array_t array = bsp_construct_array_t(dims[0], type); + bsp_error_t error = + bsp_construct_array_t_allocator(array, dims[0], type, allocator); + if (error != BSP_SUCCESS) { + H5Dclose(dset); + H5Sclose(fspace); + bsp_construct_default_array_t_allocator(array, allocator); + return BSP_ERROR_MEMORY; + } herr_t status = H5Dread(dset, bsp_get_hdf5_native_type(type), H5S_ALL, - H5S_ALL, H5P_DEFAULT, array.data); + H5S_ALL, H5P_DEFAULT, array->data); if (status < 0) { - return bsp_construct_default_array_t(); + bsp_destroy_array_t(array); + H5Sclose(fspace); + H5Dclose(dset); + bsp_construct_default_array_t_allocator(array, allocator); + return BSP_ERROR_IO; } - H5Dclose(dset); H5Sclose(fspace); - return array; + H5Dclose(dset); + return BSP_SUCCESS; } -static inline void bsp_write_attribute(hid_t f, const char* label, - const char* string) { +static inline bsp_error_t bsp_read_array(bsp_array_t* array, hid_t f, + const char* label) { + return bsp_read_array_allocator(array, f, label, bsp_default_allocator); +} + +static inline bsp_error_t bsp_write_attribute(hid_t f, const char* label, + const char* string) { hid_t strtype = H5Tcopy(H5T_C_S1); H5Tset_size(strtype, strlen(string)); H5Tset_cset(strtype, H5T_CSET_UTF8); @@ -227,23 +273,35 @@ static inline void bsp_write_attribute(hid_t f, const char* label, H5Tclose(strtype); H5Aclose(attribute); H5Sclose(dataspace); + + return BSP_SUCCESS; } -static inline char* bsp_read_attribute(hid_t f, const char* label) { +static inline bsp_error_t +bsp_read_attribute_allocator(char** string, hid_t f, const char* label, + bsp_allocator_t allocator) { hid_t attribute = H5Aopen(f, label, H5P_DEFAULT); - hid_t strtype = H5Aget_type(attribute); - hid_t type_class = H5Tget_class(strtype); - assert(type_class == H5T_STRING); + if (attribute == H5I_INVALID_HID) { + return BSP_ERROR_FORMAT; + } + + hid_t strtype = H5Aget_type(attribute); size_t size = H5Tget_size(strtype); - char* string = (char*) malloc(size + 1); + *string = (char*) allocator.malloc(size + 1); - H5Aread(attribute, strtype, string); + H5Aread(attribute, strtype, *string); + (*string)[size] = '\0'; H5Aclose(attribute); H5Tclose(strtype); - return string; + return BSP_SUCCESS; +} + +static inline bsp_error_t bsp_read_attribute(char** string, hid_t f, + const char* label) { + return bsp_read_attribute_allocator(string, f, label, bsp_default_allocator); } diff --git a/include/binsparse/matrix.h b/include/binsparse/matrix.h index db858c7..bc54f85 100644 --- a/include/binsparse/matrix.h +++ b/include/binsparse/matrix.h @@ -7,6 +7,7 @@ #pragma once #include +#include #include #include #include @@ -26,23 +27,27 @@ typedef struct bsp_matrix_t { bsp_structure_t structure; } bsp_matrix_t; -static inline bsp_matrix_t bsp_construct_default_matrix_t() { - bsp_matrix_t mat; - mat.values = bsp_construct_default_array_t(); - mat.indices_0 = bsp_construct_default_array_t(); - mat.indices_1 = bsp_construct_default_array_t(); - mat.pointers_to_1 = bsp_construct_default_array_t(); - mat.nrows = mat.ncols = mat.nnz = 0; - mat.is_iso = false; - mat.structure = BSP_GENERAL; - return mat; +static inline void +bsp_construct_default_matrix_t_allocator(bsp_matrix_t* matrix, + bsp_allocator_t allocator) { + bsp_construct_default_array_t_allocator(&matrix->values, allocator); + bsp_construct_default_array_t_allocator(&matrix->indices_0, allocator); + bsp_construct_default_array_t_allocator(&matrix->indices_1, allocator); + bsp_construct_default_array_t_allocator(&matrix->pointers_to_1, allocator); + matrix->nrows = matrix->ncols = matrix->nnz = 0; + matrix->is_iso = false; + matrix->structure = BSP_GENERAL; } -static inline void bsp_destroy_matrix_t(bsp_matrix_t matrix) { - bsp_destroy_array_t(matrix.values); - bsp_destroy_array_t(matrix.indices_0); - bsp_destroy_array_t(matrix.indices_1); - bsp_destroy_array_t(matrix.pointers_to_1); +static inline void bsp_construct_default_matrix_t(bsp_matrix_t* matrix) { + bsp_construct_default_matrix_t_allocator(matrix, bsp_default_allocator); +} + +static inline void bsp_destroy_matrix_t(bsp_matrix_t* matrix) { + bsp_destroy_array_t(&matrix->values); + bsp_destroy_array_t(&matrix->indices_0); + bsp_destroy_array_t(&matrix->indices_1); + bsp_destroy_array_t(&matrix->pointers_to_1); } static inline size_t bsp_matrix_nbytes(bsp_matrix_t mat) { diff --git a/include/binsparse/matrix_formats.h b/include/binsparse/matrix_formats.h index 2804e02..367f199 100644 --- a/include/binsparse/matrix_formats.h +++ b/include/binsparse/matrix_formats.h @@ -35,6 +35,8 @@ static inline char* bsp_get_matrix_format_string(bsp_matrix_format_t format) { return (char*) "CVEC"; } else if (format == BSP_CSR) { return (char*) "CSR"; + } else if (format == BSP_CSC) { + return (char*) "CSC"; } else if (format == BSP_DCSR) { return (char*) "DCSR"; } else if (format == BSP_DCSC) { diff --git a/include/binsparse/matrix_market/coo_sort_tools.h b/include/binsparse/matrix_market/coo_sort_tools.h index 8d3b485..a5bf934 100644 --- a/include/binsparse/matrix_market/coo_sort_tools.h +++ b/include/binsparse/matrix_market/coo_sort_tools.h @@ -45,3 +45,24 @@ static int bsp_coo_comparison_row_sort_operator_impl_(const void* x, return bsp_compare_int_impl_(x_j, y_j); } } + +static int bsp_coo_comparison_col_sort_operator_impl_(const void* x, + const void* y) { + size_t x_index = *((const size_t*) x); + size_t y_index = *((const size_t*) y); + + size_t x_i, x_j; + size_t y_i, y_j; + + bsp_array_read(bsp_coo_indices_.rowind, x_index, x_i); + bsp_array_read(bsp_coo_indices_.colind, x_index, x_j); + + bsp_array_read(bsp_coo_indices_.rowind, y_index, y_i); + bsp_array_read(bsp_coo_indices_.colind, y_index, y_j); + + if (x_j != y_j) { + return bsp_compare_int_impl_(x_j, y_j); + } else { + return bsp_compare_int_impl_(x_i, y_i); + } +} diff --git a/include/binsparse/matrix_market/matrix_market_inspector.h b/include/binsparse/matrix_market/matrix_market_inspector.h index 8811a89..0586490 100644 --- a/include/binsparse/matrix_market/matrix_market_inspector.h +++ b/include/binsparse/matrix_market/matrix_market_inspector.h @@ -8,7 +8,10 @@ #include #include +#include #include +#include +#include #include @@ -32,17 +35,38 @@ typedef struct bsp_mm_metadata { char* comments; } bsp_mm_metadata; +static inline void bsp_destroy_mm_metadata(bsp_mm_metadata* metadata) { + free(metadata->comments); +} + +static inline bool bsp_mm_metadata_is_valid(bsp_mm_metadata metadata) { + return metadata.format[0] != '\0' && metadata.type[0] != '\0' && + metadata.structure[0] != '\0'; +} + static inline bsp_mm_metadata bsp_mmread_metadata(const char* file_path) { - FILE* f = fopen(file_path, "r"); + bsp_mm_metadata metadata = {0}; + metadata.comments = (char*) calloc(1, sizeof(char)); + if (metadata.comments == NULL) { + return metadata; + } - assert(f != NULL); + FILE* f = fopen(file_path, "r"); - bsp_mm_metadata metadata; + if (f == NULL) { + return metadata; + } int read_items = fscanf(f, "%%%%MatrixMarket matrix %s %s %s\n", metadata.format, metadata.type, metadata.structure); - assert(read_items == 3); + if (read_items != 3) { + fclose(f); + metadata.format[0] = '\0'; + metadata.type[0] = '\0'; + metadata.structure[0] = '\0'; + return metadata; + } for (size_t i = 0; i < strlen(metadata.structure); i++) { metadata.structure[i] = tolower(metadata.structure[i]); @@ -54,10 +78,25 @@ static inline bsp_mm_metadata bsp_mmread_metadata(const char* file_path) { size_t comments_capacity = 2048; size_t comments_size = 0; char* comments = (char*) malloc(sizeof(char) * comments_capacity); + if (comments == NULL) { + fclose(f); + metadata.format[0] = '\0'; + metadata.type[0] = '\0'; + metadata.structure[0] = '\0'; + return metadata; + } + comments[0] = '\0'; while (!outOfComments) { char* line = fgets(buf, 2048, f); - assert(line != NULL); + if (line == NULL) { + free(comments); + fclose(f); + metadata.format[0] = '\0'; + metadata.type[0] = '\0'; + metadata.structure[0] = '\0'; + return metadata; + } if (line[0] != '%') { outOfComments = 1; @@ -68,25 +107,49 @@ static inline bsp_mm_metadata bsp_mmread_metadata(const char* file_path) { while (comments_size + strlen(line) > comments_capacity) { comments_capacity <<= 1; } - comments = (char*) realloc(comments, sizeof(char) * comments_capacity); + char* resized = + (char*) realloc(comments, sizeof(char) * comments_capacity); + if (resized == NULL) { + free(comments); + fclose(f); + metadata.format[0] = '\0'; + metadata.type[0] = '\0'; + metadata.structure[0] = '\0'; + return metadata; + } + comments = resized; } memcpy(comments + comments_size, line, strlen(line)); comments_size += strlen(line); + comments[comments_size] = '\0'; } } - if (comments[comments_size - 1] == '\n') { + if (comments_size > 0 && comments[comments_size - 1] == '\n') { comments[comments_size - 1] = 0; } + free(metadata.comments); metadata.comments = comments; unsigned long long nrows, ncols, nnz; if (strcmp(metadata.format, "coordinate") == 0) { - sscanf(buf, "%llu %llu %llu", &nrows, &ncols, &nnz); + if (sscanf(buf, "%llu %llu %llu", &nrows, &ncols, &nnz) != 3) { + fclose(f); + metadata.format[0] = '\0'; + metadata.type[0] = '\0'; + metadata.structure[0] = '\0'; + return metadata; + } } else { - sscanf(buf, "%llu %llu", &nrows, &ncols); + if (sscanf(buf, "%llu %llu", &nrows, &ncols) != 2) { + fclose(f); + metadata.format[0] = '\0'; + metadata.type[0] = '\0'; + metadata.structure[0] = '\0'; + return metadata; + } nnz = nrows * ncols; } diff --git a/include/binsparse/matrix_market/matrix_market_read.h b/include/binsparse/matrix_market/matrix_market_read.h index c823a40..39971d3 100644 --- a/include/binsparse/matrix_market/matrix_market_read.h +++ b/include/binsparse/matrix_market/matrix_market_read.h @@ -32,13 +32,18 @@ static inline bsp_matrix_t bsp_mmread_explicit_array(const char* file_path, assert(false); } - bsp_matrix_t matrix = bsp_construct_default_matrix_t(); + bsp_matrix_t matrix; + bsp_construct_default_matrix_t(&matrix); matrix.nrows = metadata.nrows; matrix.ncols = metadata.ncols; matrix.nnz = matrix.nrows * matrix.ncols; - matrix.values = bsp_construct_array_t(matrix.nnz, value_type); + bsp_error_t error = + bsp_construct_array_t(&matrix.values, matrix.nnz, value_type); + if (error != BSP_SUCCESS) { + return matrix; // Return empty matrix on error + } matrix.format = BSP_DMAT; @@ -104,6 +109,8 @@ static inline bsp_matrix_t bsp_mmread_explicit_array(const char* file_path, fclose(f); + bsp_destroy_mm_metadata(&metadata); + return matrix; } @@ -125,7 +132,8 @@ bsp_mmread_explicit_coordinate(const char* file_path, bsp_type_t value_type, assert(false); } - bsp_matrix_t matrix = bsp_construct_default_matrix_t(); + bsp_matrix_t matrix; + bsp_construct_default_matrix_t(&matrix); if (mm_type == BSP_MM_PATTERN) { matrix.is_iso = true; @@ -135,14 +143,33 @@ bsp_mmread_explicit_coordinate(const char* file_path, bsp_type_t value_type, matrix.ncols = metadata.ncols; matrix.nnz = metadata.nnz; - matrix.indices_0 = bsp_construct_array_t(matrix.nnz, index_type); - matrix.indices_1 = bsp_construct_array_t(matrix.nnz, index_type); + bsp_error_t error = + bsp_construct_array_t(&matrix.indices_0, matrix.nnz, index_type); + if (error != BSP_SUCCESS) { + return matrix; + } + + error = bsp_construct_array_t(&matrix.indices_1, matrix.nnz, index_type); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&matrix.indices_0); + return matrix; + } if (mm_type == BSP_MM_PATTERN) { - matrix.values = bsp_construct_array_t(1, value_type); + error = bsp_construct_array_t(&matrix.values, 1, value_type); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&matrix.indices_0); + bsp_destroy_array_t(&matrix.indices_1); + return matrix; + } bsp_array_write(matrix.values, 0, true); } else { - matrix.values = bsp_construct_array_t(matrix.nnz, value_type); + error = bsp_construct_array_t(&matrix.values, matrix.nnz, value_type); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&matrix.indices_0); + bsp_destroy_array_t(&matrix.indices_1); + return matrix; + } } matrix.format = BSP_COO; @@ -235,13 +262,29 @@ bsp_mmread_explicit_coordinate(const char* file_path, bsp_type_t value_type, qsort(indices, matrix.nnz, sizeof(size_t), bsp_coo_comparison_row_sort_operator_impl_); - bsp_array_t rowind = bsp_copy_construct_array_t(matrix.indices_0); - bsp_array_t colind = bsp_copy_construct_array_t(matrix.indices_1); + bsp_array_t rowind; + bsp_array_t colind; + + error = bsp_copy_construct_array_t(&rowind, matrix.indices_0); + if (error != BSP_SUCCESS) { + return matrix; + } + + error = bsp_copy_construct_array_t(&colind, matrix.indices_1); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&rowind); + return matrix; + } bsp_array_t values; if (!matrix.is_iso) { - values = bsp_copy_construct_array_t(matrix.values); + error = bsp_copy_construct_array_t(&values, matrix.values); + if (error != BSP_SUCCESS) { + bsp_destroy_array_t(&rowind); + bsp_destroy_array_t(&colind); + return matrix; + } } for (size_t i = 0; i < matrix.nnz; i++) { @@ -252,19 +295,21 @@ bsp_mmread_explicit_coordinate(const char* file_path, bsp_type_t value_type, } } - bsp_destroy_array_t(matrix.indices_0); - bsp_destroy_array_t(matrix.indices_1); + bsp_destroy_array_t(&matrix.indices_0); + bsp_destroy_array_t(&matrix.indices_1); matrix.indices_0 = rowind; matrix.indices_1 = colind; if (!matrix.is_iso) { - bsp_destroy_array_t(matrix.values); + bsp_destroy_array_t(&matrix.values); matrix.values = values; } free(indices); fclose(f); + bsp_destroy_mm_metadata(&metadata); + return matrix; } @@ -272,10 +317,18 @@ static inline bsp_matrix_t bsp_mmread_explicit(const char* file_path, bsp_type_t value_type, bsp_type_t index_type) { bsp_mm_metadata metadata = bsp_mmread_metadata(file_path); + if (!bsp_mm_metadata_is_valid(metadata)) { + bsp_matrix_t matrix; + bsp_construct_default_matrix_t(&matrix); + bsp_destroy_mm_metadata(&metadata); + return matrix; + } if (strcmp(metadata.format, "array") == 0) { + bsp_destroy_mm_metadata(&metadata); return bsp_mmread_explicit_array(file_path, value_type, index_type); } else if (strcmp(metadata.format, "coordinate") == 0) { + bsp_destroy_mm_metadata(&metadata); return bsp_mmread_explicit_coordinate(file_path, value_type, index_type); } else { assert(false); @@ -284,6 +337,12 @@ static inline bsp_matrix_t bsp_mmread_explicit(const char* file_path, static inline bsp_matrix_t bsp_mmread(const char* file_path) { bsp_mm_metadata metadata = bsp_mmread_metadata(file_path); + if (!bsp_mm_metadata_is_valid(metadata)) { + bsp_matrix_t matrix; + bsp_construct_default_matrix_t(&matrix); + bsp_destroy_mm_metadata(&metadata); + return matrix; + } bsp_type_t value_type; @@ -304,5 +363,7 @@ static inline bsp_matrix_t bsp_mmread(const char* file_path) { bsp_type_t index_type = bsp_pick_integer_type(max_dim); + bsp_destroy_mm_metadata(&metadata); + return bsp_mmread_explicit(file_path, value_type, index_type); } diff --git a/include/binsparse/minimize_values.h b/include/binsparse/minimize_values.h index 4969000..7ef18f8 100644 --- a/include/binsparse/minimize_values.h +++ b/include/binsparse/minimize_values.h @@ -24,8 +24,12 @@ static inline bsp_matrix_t bsp_matrix_minimize_values(bsp_matrix_t matrix) { } if (float32_representable) { - bsp_array_t new_values = - bsp_construct_array_t(matrix.values.size, BSP_FLOAT32); + bsp_array_t new_values; + bsp_error_t error = + bsp_construct_array_t(&new_values, matrix.values.size, BSP_FLOAT32); + if (error != BSP_SUCCESS) { + return matrix; // Return original matrix on error + } float* n_values = (float*) new_values.data; @@ -33,7 +37,7 @@ static inline bsp_matrix_t bsp_matrix_minimize_values(bsp_matrix_t matrix) { n_values[i] = values[i]; } - bsp_destroy_array_t(matrix.values); + bsp_destroy_array_t(&matrix.values); matrix.values = new_values; } } else if (matrix.values.type == BSP_INT64) { @@ -78,8 +82,12 @@ static inline bsp_matrix_t bsp_matrix_minimize_values(bsp_matrix_t matrix) { value_type = BSP_INT64; } } - bsp_array_t new_values = - bsp_construct_array_t(matrix.values.size, value_type); + bsp_array_t new_values; + bsp_error_t error = + bsp_construct_array_t(&new_values, matrix.values.size, value_type); + if (error != BSP_SUCCESS) { + return matrix; // Return original matrix on error + } for (size_t i = 0; i < matrix.values.size; i++) { int64_t value; @@ -87,7 +95,7 @@ static inline bsp_matrix_t bsp_matrix_minimize_values(bsp_matrix_t matrix) { bsp_array_write(new_values, i, value); } - bsp_destroy_array_t(matrix.values); + bsp_destroy_array_t(&matrix.values); matrix.values = new_values; } else if (matrix.values.type == BSP_COMPLEX_FLOAT64) { bool float32_representable = true; @@ -101,8 +109,12 @@ static inline bsp_matrix_t bsp_matrix_minimize_values(bsp_matrix_t matrix) { } if (float32_representable) { - bsp_array_t new_values = - bsp_construct_array_t(matrix.values.size, BSP_COMPLEX_FLOAT32); + bsp_array_t new_values; + bsp_error_t error = bsp_construct_array_t(&new_values, matrix.values.size, + BSP_COMPLEX_FLOAT32); + if (error != BSP_SUCCESS) { + return matrix; // Return original matrix on error + } float _Complex* n_values = (float _Complex*) new_values.data; @@ -110,7 +122,7 @@ static inline bsp_matrix_t bsp_matrix_minimize_values(bsp_matrix_t matrix) { n_values[i] = values[i]; } - bsp_destroy_array_t(matrix.values); + bsp_destroy_array_t(&matrix.values); matrix.values = new_values; } } diff --git a/include/binsparse/read_matrix.h b/include/binsparse/read_matrix.h index 948ffba..523f34c 100644 --- a/include/binsparse/read_matrix.h +++ b/include/binsparse/read_matrix.h @@ -6,26 +6,24 @@ #pragma once +#include +#include + #ifdef __cplusplus extern "C" { #endif -#ifdef BSP_USE_HDF5 -#include - -#if __STDC_VERSION__ >= 201112L -bsp_matrix_t bsp_read_matrix_from_group_parallel(hid_t f, int num_threads); -#endif - -bsp_matrix_t bsp_read_matrix_from_group(hid_t f); -#endif - #if __STDC_VERSION__ >= 201112L -bsp_matrix_t bsp_read_matrix_parallel(const char* file_name, const char* group, - int num_threads); +bsp_error_t bsp_read_matrix_parallel(bsp_matrix_t* matrix, + const char* file_name, const char* group, + int num_threads); #endif -bsp_matrix_t bsp_read_matrix(const char* file_name, const char* group); +bsp_error_t bsp_read_matrix(bsp_matrix_t* matrix, const char* file_name, + const char* group); +bsp_error_t bsp_read_matrix_allocator(bsp_matrix_t* matrix, + const char* file_name, const char* group, + bsp_allocator_t allocator); #ifdef __cplusplus } diff --git a/include/binsparse/tensor.h b/include/binsparse/tensor.h index 4d63512..31d6b89 100644 --- a/include/binsparse/tensor.h +++ b/include/binsparse/tensor.h @@ -74,7 +74,7 @@ static void bsp_destroy_level_t(bsp_level_t* level) { switch (level->kind) { case BSP_TENSOR_ELEMENT: { bsp_element_t* element = (bsp_element_t*) level->data; - bsp_destroy_array_t(element->values); + bsp_destroy_array_t(&element->values); free(element); break; } @@ -88,10 +88,10 @@ static void bsp_destroy_level_t(bsp_level_t* level) { bsp_sparse_t* sparse = (bsp_sparse_t*) level->data; if (sparse->pointers_to != NULL) - bsp_destroy_array_t(*sparse->pointers_to); + bsp_destroy_array_t(sparse->pointers_to); if (sparse->indices != NULL) { for (int i = 0; i < sparse->rank; i++) { - bsp_destroy_array_t(sparse->indices[i]); + bsp_destroy_array_t(&sparse->indices[i]); } } bsp_destroy_level_t(sparse->child); diff --git a/include/binsparse/types.h b/include/binsparse/types.h index e3c9e53..a98e938 100644 --- a/include/binsparse/types.h +++ b/include/binsparse/types.h @@ -7,7 +7,8 @@ #pragma once #include -#include +#include +#include typedef enum bsp_type_t { BSP_UINT8 = 0, @@ -90,212 +91,14 @@ static inline size_t bsp_type_size(bsp_type_t type) { } } -static inline hid_t bsp_get_hdf5_standard_type(bsp_type_t type) { - if (type == BSP_UINT8) { - return H5T_STD_U8LE; - } else if (type == BSP_UINT16) { - return H5T_STD_U16LE; - } else if (type == BSP_UINT32) { - return H5T_STD_U32LE; - } else if (type == BSP_UINT64) { - return H5T_STD_U64LE; - } else if (type == BSP_INT8) { - return H5T_STD_I8LE; - } else if (type == BSP_INT16) { - return H5T_STD_I16LE; - } else if (type == BSP_INT32) { - return H5T_STD_I32LE; - } else if (type == BSP_INT64) { - return H5T_STD_I64LE; - } else if (type == BSP_FLOAT32) { - return H5T_IEEE_F32LE; - } else if (type == BSP_FLOAT64) { - return H5T_IEEE_F64LE; - } else if (type == BSP_BINT8) { - return H5T_STD_I8LE; - } else { - return H5I_INVALID_HID; - } -} - -static inline bsp_type_t bsp_get_bsp_type(hid_t type) { - H5T_class_t cl = H5Tget_class(type); - H5T_order_t order = H5Tget_order(type); - H5T_sign_t sign = H5Tget_sign(type); - size_t size = H5Tget_size(type); - - if (cl == H5T_INTEGER) { - if (sign == H5T_SGN_NONE) { - if (size == 1) { - return BSP_UINT8; - } else if (size == 2) { - return BSP_UINT16; - } else if (size == 4) { - return BSP_UINT32; - } else if (size == 8) { - return BSP_UINT64; - } else { - return BSP_INVALID_TYPE; - } - } else /* if (sign == H5T_SGN_2) */ { - if (size == 1) { - return BSP_INT8; - } else if (size == 2) { - return BSP_INT16; - } else if (size == 4) { - return BSP_INT32; - } else if (size == 8) { - return BSP_INT64; - } else { - return BSP_INVALID_TYPE; - } - } - } else if (cl == H5T_FLOAT) { - if (size == 4) { - return BSP_FLOAT32; - } else if (size == 8) { - return BSP_FLOAT64; - } else { - return BSP_INVALID_TYPE; - } - } else { - return BSP_INVALID_TYPE; - } -} - -// NOTE: This code is a bit silly, but it seems to be the only -// way to generically determine the HDF5 native types for -// stdint's fixed width integer types. -static inline hid_t bsp_get_hdf5_native_type(bsp_type_t type) { - if (type == BSP_INT8 || type == BSP_BINT8) { - if (sizeof(int8_t) == sizeof(char)) { - return H5T_NATIVE_CHAR; - } else if (sizeof(int8_t) == sizeof(short)) { - return H5T_NATIVE_SHORT; - } else if (sizeof(int8_t) == sizeof(int)) { - return H5T_NATIVE_INT; - } else if (sizeof(int8_t) == sizeof(long)) { - return H5T_NATIVE_LONG; - } else if (sizeof(int8_t) == sizeof(long long)) { - return H5T_NATIVE_LLONG; - } else { - assert(false); - } - } else if (type == BSP_INT16) { - if (sizeof(int16_t) == sizeof(char)) { - return H5T_NATIVE_CHAR; - } else if (sizeof(int16_t) == sizeof(short)) { - return H5T_NATIVE_SHORT; - } else if (sizeof(int16_t) == sizeof(int)) { - return H5T_NATIVE_INT; - } else if (sizeof(int32_t) == sizeof(long)) { - return H5T_NATIVE_LONG; - } else if (sizeof(int64_t) == sizeof(long long)) { - return H5T_NATIVE_LLONG; - } else { - assert(false); - } - } else if (type == BSP_INT32) { - if (sizeof(int32_t) == sizeof(char)) { - return H5T_NATIVE_CHAR; - } else if (sizeof(int32_t) == sizeof(short)) { - return H5T_NATIVE_SHORT; - } else if (sizeof(int32_t) == sizeof(int)) { - return H5T_NATIVE_INT; - } else if (sizeof(int32_t) == sizeof(long)) { - return H5T_NATIVE_LONG; - } else if (sizeof(int32_t) == sizeof(long long)) { - return H5T_NATIVE_LLONG; - } else { - assert(false); - } - } else if (type == BSP_INT64) { - if (sizeof(int64_t) == sizeof(char)) { - return H5T_NATIVE_CHAR; - } else if (sizeof(int64_t) == sizeof(short)) { - return H5T_NATIVE_SHORT; - } else if (sizeof(int64_t) == sizeof(int)) { - return H5T_NATIVE_INT; - } else if (sizeof(int64_t) == sizeof(long)) { - return H5T_NATIVE_LONG; - } else if (sizeof(int64_t) == sizeof(long long)) { - return H5T_NATIVE_LLONG; - } else { - assert(false); - } - } else if (type == BSP_UINT8) { - if (sizeof(uint8_t) == sizeof(unsigned char)) { - return H5T_NATIVE_UCHAR; - } else if (sizeof(uint8_t) == sizeof(unsigned short)) { - return H5T_NATIVE_USHORT; - } else if (sizeof(uint8_t) == sizeof(unsigned int)) { - return H5T_NATIVE_UINT; - } else if (sizeof(uint8_t) == sizeof(unsigned long)) { - return H5T_NATIVE_ULONG; - } else if (sizeof(uint8_t) == sizeof(unsigned long long)) { - return H5T_NATIVE_ULLONG; - } else { - assert(false); - } - } else if (type == BSP_UINT16) { - if (sizeof(uint16_t) == sizeof(unsigned char)) { - return H5T_NATIVE_UCHAR; - } else if (sizeof(uint16_t) == sizeof(unsigned short)) { - return H5T_NATIVE_USHORT; - } else if (sizeof(uint16_t) == sizeof(unsigned int)) { - return H5T_NATIVE_UINT; - } else if (sizeof(uint16_t) == sizeof(unsigned long)) { - return H5T_NATIVE_ULONG; - } else if (sizeof(uint16_t) == sizeof(unsigned long long)) { - return H5T_NATIVE_ULLONG; - } else { - assert(false); - } - } else if (type == BSP_UINT32) { - if (sizeof(uint32_t) == sizeof(unsigned char)) { - return H5T_NATIVE_UCHAR; - } else if (sizeof(uint32_t) == sizeof(unsigned short)) { - return H5T_NATIVE_USHORT; - } else if (sizeof(uint32_t) == sizeof(unsigned int)) { - return H5T_NATIVE_UINT; - } else if (sizeof(uint32_t) == sizeof(unsigned long)) { - return H5T_NATIVE_ULONG; - } else if (sizeof(uint32_t) == sizeof(unsigned long long)) { - return H5T_NATIVE_ULLONG; - } else { - assert(false); - } - } else if (type == BSP_UINT64) { - if (sizeof(uint64_t) == sizeof(unsigned char)) { - return H5T_NATIVE_UCHAR; - } else if (sizeof(uint64_t) == sizeof(unsigned short)) { - return H5T_NATIVE_USHORT; - } else if (sizeof(uint64_t) == sizeof(unsigned int)) { - return H5T_NATIVE_UINT; - } else if (sizeof(uint64_t) == sizeof(unsigned long)) { - return H5T_NATIVE_ULONG; - } else if (sizeof(uint64_t) == sizeof(unsigned long long)) { - return H5T_NATIVE_ULLONG; - } else { - assert(false); - } - } else if (type == BSP_FLOAT32) { - return H5T_NATIVE_FLOAT; - } else if (type == BSP_FLOAT64) { - return H5T_NATIVE_DOUBLE; - } else { - return H5I_INVALID_HID; - } -} - // Given the maximum value `max_value` that must be stored, // pick an unsigned integer type for indices. static inline bsp_type_t bsp_pick_integer_type(size_t max_value) { - if (max_value < (size_t) UINT8_MAX) { + if (max_value <= (size_t) UINT8_MAX) { return BSP_UINT8; - } else if (max_value < (size_t) UINT16_MAX) { + } else if (max_value <= (size_t) UINT16_MAX) { return BSP_UINT16; - } else if (max_value < (size_t) UINT32_MAX) { + } else if (max_value <= (size_t) UINT32_MAX) { return BSP_UINT32; } else { return BSP_UINT64; diff --git a/include/binsparse/write_matrix.h b/include/binsparse/write_matrix.h index 53b89f5..e535569 100644 --- a/include/binsparse/write_matrix.h +++ b/include/binsparse/write_matrix.h @@ -13,17 +13,10 @@ extern "C" { // TODO: make cJSON optional. #include -#include -#ifdef BSP_USE_HDF5 -#include - -int bsp_write_matrix_to_group(hid_t f, bsp_matrix_t matrix, cJSON* user_json, - int compression_level); -#endif - -int bsp_write_matrix(const char* fname, bsp_matrix_t matrix, const char* group, - cJSON* user_json, int compression_level); +bsp_error_t bsp_write_matrix(const char* fname, bsp_matrix_t matrix, + const char* group, const char* user_json, + int compression_level); #ifdef __cplusplus } diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 0b82bd2..f966494 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -8,3 +8,10 @@ target_sources(binsparse PRIVATE write_matrix.c write_tensor.c ) + +target_sources(binsparse_dynamic PRIVATE + read_matrix.c + read_tensor.c + write_matrix.c + write_tensor.c +) diff --git a/src/read_matrix.c b/src/read_matrix.c index 464210c..103851b 100644 --- a/src/read_matrix.c +++ b/src/read_matrix.c @@ -5,21 +5,43 @@ */ #include +#include +#include +#include #include #include #include #include +#include + +static void bsp_prepare_hdf5_runtime(void) { + static bool initialized = false; + if (!initialized) { + H5dont_atexit(); + initialized = true; + } +} #if __STDC_VERSION__ >= 201112L -bsp_matrix_t bsp_read_matrix_from_group_parallel(hid_t f, int num_threads) { - bsp_matrix_t matrix = bsp_construct_default_matrix_t(); +bsp_error_t bsp_read_matrix_from_group_parallel(bsp_matrix_t* matrix, hid_t f, + int num_threads) { + if (f < 0) { + return BSP_ERROR_IO; + } - char* json_string = bsp_read_attribute(f, (char*) "binsparse"); + bsp_construct_default_matrix_t(matrix); + + char* json_string; + bsp_error_t error = bsp_read_attribute(&json_string, f, (char*) "binsparse"); + if (error != BSP_SUCCESS) { + return error; + } cJSON* j = cJSON_Parse(json_string); - assert(j != NULL); - assert(cJSON_IsObject(j)); + if (j == NULL || !cJSON_IsObject(j)) { + return BSP_ERROR_FORMAT; + } cJSON* binsparse = cJSON_GetObjectItemCaseSensitive(j, "binsparse"); assert(cJSON_IsObject(binsparse)); @@ -40,7 +62,7 @@ bsp_matrix_t bsp_read_matrix_from_group_parallel(hid_t f, int num_threads) { assert(format != 0); - matrix.format = format; + matrix->format = format; cJSON* nnz_ = cJSON_GetObjectItemCaseSensitive(binsparse, "number_of_stored_values"); @@ -62,159 +84,265 @@ bsp_matrix_t bsp_read_matrix_from_group_parallel(hid_t f, int num_threads) { size_t ncols = cJSON_GetNumberValue(ncols_); - matrix.nrows = nrows; - matrix.ncols = ncols; - matrix.nnz = nnz; - matrix.format = format; + matrix->nrows = nrows; + matrix->ncols = ncols; + matrix->nnz = nnz; + matrix->format = format; cJSON* data_types_ = cJSON_GetObjectItemCaseSensitive(binsparse, "data_types"); assert(data_types_ != NULL); if (cJSON_HasObjectItem(data_types_, "values")) { - matrix.values = bsp_read_array_parallel(f, (char*) "values", num_threads); + error = bsp_read_array_parallel(&matrix->values, f, (char*) "values", + num_threads); + if (error != BSP_SUCCESS) { + free(json_string); + return error; + } cJSON* value_type = cJSON_GetObjectItemCaseSensitive(data_types_, "values"); char* type_string = cJSON_GetStringValue(value_type); if (strlen(type_string) >= 4 && strncmp(type_string, "iso[", 4) == 0) { - matrix.is_iso = true; + matrix->is_iso = true; type_string += 4; } if (strlen(type_string) >= 8 && strncmp(type_string, "complex[", 8) == 0) { - matrix.values = bsp_fp_array_to_complex(matrix.values); + bsp_error_t error = bsp_fp_array_to_complex(&matrix->values); + if (error != BSP_SUCCESS) { + // TODO: handle error + return error; + } } } if (cJSON_HasObjectItem(data_types_, "indices_0")) { - matrix.indices_0 = - bsp_read_array_parallel(f, (char*) "indices_0", num_threads); + error = bsp_read_array_parallel(&matrix->indices_0, f, (char*) "indices_0", + num_threads); + if (error != BSP_SUCCESS) { + free(json_string); + bsp_destroy_array_t(&matrix->values); + return error; + } } if (cJSON_HasObjectItem(data_types_, "indices_1")) { - matrix.indices_1 = - bsp_read_array_parallel(f, (char*) "indices_1", num_threads); + error = bsp_read_array_parallel(&matrix->indices_1, f, (char*) "indices_1", + num_threads); + if (error != BSP_SUCCESS) { + free(json_string); + bsp_destroy_array_t(&matrix->values); + bsp_destroy_array_t(&matrix->indices_0); + return error; + } } if (cJSON_HasObjectItem(data_types_, "pointers_to_1")) { - matrix.pointers_to_1 = - bsp_read_array_parallel(f, (char*) "pointers_to_1", num_threads); + error = bsp_read_array_parallel(&matrix->pointers_to_1, f, + (char*) "pointers_to_1", num_threads); + if (error != BSP_SUCCESS) { + free(json_string); + bsp_destroy_array_t(&matrix->values); + bsp_destroy_array_t(&matrix->indices_0); + bsp_destroy_array_t(&matrix->indices_1); + return error; + } } if (cJSON_HasObjectItem(binsparse, "structure")) { cJSON* structure_ = cJSON_GetObjectItemCaseSensitive(binsparse, "structure"); char* structure = cJSON_GetStringValue(structure_); - matrix.structure = bsp_get_structure(structure); + matrix->structure = bsp_get_structure(structure); } cJSON_Delete(j); free(json_string); - return matrix; + return BSP_SUCCESS; } #endif -bsp_matrix_t bsp_read_matrix_from_group(hid_t f) { - bsp_matrix_t matrix = bsp_construct_default_matrix_t(); +bsp_error_t bsp_read_matrix_from_group_allocator(bsp_matrix_t* matrix, hid_t f, + bsp_allocator_t allocator) { + if (f < 0) { + return BSP_ERROR_IO; + } + + bsp_construct_default_matrix_t_allocator(matrix, allocator); - char* json_string = bsp_read_attribute(f, (char*) "binsparse"); + char* json_string; + bsp_error_t error = bsp_read_attribute_allocator( + &json_string, f, (char*) "binsparse", allocator); + if (error != BSP_SUCCESS) { + return error; + } cJSON* j = cJSON_Parse(json_string); - assert(j != NULL); - assert(cJSON_IsObject(j)); + if (j == NULL || !cJSON_IsObject(j)) { + allocator.free(json_string); + return BSP_ERROR_FORMAT; + } cJSON* binsparse = cJSON_GetObjectItemCaseSensitive(j, "binsparse"); - assert(cJSON_IsObject(binsparse)); + if (!cJSON_IsObject(binsparse)) { + cJSON_Delete(j); + allocator.free(json_string); + return BSP_ERROR_FORMAT; + } cJSON* version_ = cJSON_GetObjectItemCaseSensitive(binsparse, "version"); - assert(version_ != NULL); - - assert(cJSON_IsString(version_)); + if (version_ == NULL || !cJSON_IsString(version_)) { + cJSON_Delete(j); + allocator.free(json_string); + return BSP_ERROR_FORMAT; + } // TODO: check version. cJSON* format_ = cJSON_GetObjectItemCaseSensitive(binsparse, "format"); - assert(format_ != NULL); + if (format_ == NULL || !cJSON_IsString(format_)) { + cJSON_Delete(j); + allocator.free(json_string); + return BSP_ERROR_FORMAT; + } char* format_string = cJSON_GetStringValue(format_); bsp_matrix_format_t format = bsp_get_matrix_format(format_string); - assert(format != 0); + if (format == BSP_INVALID_FORMAT) { + cJSON_Delete(j); + allocator.free(json_string); + return BSP_ERROR_FORMAT; + } - matrix.format = format; + matrix->format = format; cJSON* nnz_ = cJSON_GetObjectItemCaseSensitive(binsparse, "number_of_stored_values"); - assert(nnz_ != NULL); + if (nnz_ == NULL || !cJSON_IsNumber(nnz_)) { + cJSON_Delete(j); + allocator.free(json_string); + return BSP_ERROR_FORMAT; + } size_t nnz = cJSON_GetNumberValue(nnz_); cJSON* shape_ = cJSON_GetObjectItemCaseSensitive(binsparse, "shape"); - assert(shape_ != NULL); - - assert(cJSON_GetArraySize(shape_) == 2); + if (shape_ == NULL || !cJSON_IsArray(shape_) || + cJSON_GetArraySize(shape_) != 2) { + cJSON_Delete(j); + allocator.free(json_string); + return BSP_ERROR_FORMAT; + } cJSON* nrows_ = cJSON_GetArrayItem(shape_, 0); - assert(nrows_ != NULL); + if (nrows_ == NULL || !cJSON_IsNumber(nrows_)) { + cJSON_Delete(j); + allocator.free(json_string); + return BSP_ERROR_FORMAT; + } size_t nrows = cJSON_GetNumberValue(nrows_); cJSON* ncols_ = cJSON_GetArrayItem(shape_, 1); - assert(ncols_ != NULL); + if (ncols_ == NULL || !cJSON_IsNumber(ncols_)) { + cJSON_Delete(j); + allocator.free(json_string); + return BSP_ERROR_FORMAT; + } size_t ncols = cJSON_GetNumberValue(ncols_); - matrix.nrows = nrows; - matrix.ncols = ncols; - matrix.nnz = nnz; - matrix.format = format; + matrix->nrows = nrows; + matrix->ncols = ncols; + matrix->nnz = nnz; + matrix->format = format; cJSON* data_types_ = cJSON_GetObjectItemCaseSensitive(binsparse, "data_types"); - assert(data_types_ != NULL); + if (data_types_ == NULL || !cJSON_IsObject(data_types_)) { + cJSON_Delete(j); + allocator.free(json_string); + return BSP_ERROR_FORMAT; + } if (cJSON_HasObjectItem(data_types_, "values")) { - matrix.values = bsp_read_array(f, (char*) "values"); + error = bsp_read_array_allocator(&matrix->values, f, (char*) "values", + allocator); + if (error != BSP_SUCCESS) { + allocator.free(json_string); + return error; + } cJSON* value_type = cJSON_GetObjectItemCaseSensitive(data_types_, "values"); char* type_string = cJSON_GetStringValue(value_type); if (strlen(type_string) >= 4 && strncmp(type_string, "iso[", 4) == 0) { - matrix.is_iso = true; + matrix->is_iso = true; type_string += 4; } if (strlen(type_string) >= 8 && strncmp(type_string, "complex[", 8) == 0) { - matrix.values = bsp_fp_array_to_complex(matrix.values); + bsp_error_t error = bsp_fp_array_to_complex(&matrix->values); + if (error != BSP_SUCCESS) { + // TODO: handle error + return error; + } } } if (cJSON_HasObjectItem(data_types_, "indices_0")) { - matrix.indices_0 = bsp_read_array(f, (char*) "indices_0"); + error = bsp_read_array_allocator(&matrix->indices_0, f, (char*) "indices_0", + allocator); + if (error != BSP_SUCCESS) { + allocator.free(json_string); + bsp_destroy_array_t(&matrix->values); + return error; + } } if (cJSON_HasObjectItem(data_types_, "indices_1")) { - matrix.indices_1 = bsp_read_array(f, (char*) "indices_1"); + error = bsp_read_array_allocator(&matrix->indices_1, f, (char*) "indices_1", + allocator); + if (error != BSP_SUCCESS) { + allocator.free(json_string); + bsp_destroy_array_t(&matrix->values); + bsp_destroy_array_t(&matrix->indices_0); + return error; + } } if (cJSON_HasObjectItem(data_types_, "pointers_to_1")) { - matrix.pointers_to_1 = bsp_read_array(f, (char*) "pointers_to_1"); + error = bsp_read_array_allocator(&matrix->pointers_to_1, f, + (char*) "pointers_to_1", allocator); + if (error != BSP_SUCCESS) { + allocator.free(json_string); + bsp_destroy_array_t(&matrix->values); + bsp_destroy_array_t(&matrix->indices_0); + bsp_destroy_array_t(&matrix->indices_1); + return error; + } } if (cJSON_HasObjectItem(binsparse, "structure")) { cJSON* structure_ = cJSON_GetObjectItemCaseSensitive(binsparse, "structure"); char* structure = cJSON_GetStringValue(structure_); - matrix.structure = bsp_get_structure(structure); + matrix->structure = bsp_get_structure(structure); } cJSON_Delete(j); - free(json_string); + allocator.free(json_string); + + return BSP_SUCCESS; +} - return matrix; +bsp_error_t bsp_read_matrix_from_group(bsp_matrix_t* matrix, hid_t f) { + return bsp_read_matrix_from_group_allocator(matrix, f, bsp_default_allocator); } static inline size_t bsp_final_dot(const char* str) { @@ -228,52 +356,105 @@ static inline size_t bsp_final_dot(const char* str) { } #if __STDC_VERSION__ >= 201112L -bsp_matrix_t bsp_read_matrix_parallel(const char* file_name, const char* group, - int num_threads) { +bsp_error_t bsp_read_matrix_parallel(bsp_matrix_t* matrix, + const char* file_name, const char* group, + int num_threads) { if (group == NULL) { size_t idx = bsp_final_dot(file_name); if (strcmp(file_name + idx, ".hdf5") == 0 || strcmp(file_name + idx, ".h5") == 0) { + if (access(file_name, F_OK) != 0) { + return BSP_ERROR_IO; + } + bsp_prepare_hdf5_runtime(); hid_t f = H5Fopen(file_name, H5F_ACC_RDONLY, H5P_DEFAULT); - bsp_matrix_t matrix = bsp_read_matrix_from_group_parallel(f, num_threads); + if (f < 0) { + return BSP_ERROR_IO; + } + bsp_error_t error = + bsp_read_matrix_from_group_parallel(matrix, f, num_threads); H5Fclose(f); - return matrix; + return error; } else if (strcmp(file_name + idx, ".mtx") == 0) { - return bsp_mmread(file_name); + // TODO: implement error-handling for Matrix Market. + *matrix = bsp_mmread(file_name); + return BSP_SUCCESS; } else { - assert(false); + return BSP_ERROR_IO; } } else { + if (access(file_name, F_OK) != 0) { + return BSP_ERROR_IO; + } + bsp_prepare_hdf5_runtime(); hid_t f = H5Fopen(file_name, H5F_ACC_RDONLY, H5P_DEFAULT); + if (f < 0) { + return BSP_ERROR_IO; + } hid_t g = H5Gopen1(f, group); - bsp_matrix_t matrix = bsp_read_matrix_from_group_parallel(g, num_threads); + if (g < 0) { + H5Fclose(f); + return BSP_ERROR_IO; + } + bsp_error_t error = + bsp_read_matrix_from_group_parallel(matrix, g, num_threads); H5Gclose(g); H5Fclose(f); - return matrix; + return error; } } #endif -bsp_matrix_t bsp_read_matrix(const char* file_name, const char* group) { +bsp_error_t bsp_read_matrix_allocator(bsp_matrix_t* matrix, + const char* file_name, const char* group, + bsp_allocator_t allocator) { if (group == NULL) { size_t idx = bsp_final_dot(file_name); if (strcmp(file_name + idx, ".hdf5") == 0 || strcmp(file_name + idx, ".h5") == 0) { + if (access(file_name, F_OK) != 0) { + return BSP_ERROR_IO; + } + bsp_prepare_hdf5_runtime(); hid_t f = H5Fopen(file_name, H5F_ACC_RDONLY, H5P_DEFAULT); - bsp_matrix_t matrix = bsp_read_matrix_from_group(f); + if (f < 0) { + return BSP_ERROR_IO; + } + bsp_error_t error = + bsp_read_matrix_from_group_allocator(matrix, f, allocator); H5Fclose(f); - return matrix; + return error; } else if (strcmp(file_name + idx, ".mtx") == 0) { - return bsp_mmread(file_name); + // TODO: implement error-handling for Matrix Market. + *matrix = bsp_mmread(file_name); + return BSP_SUCCESS; } else { - assert(false); + return BSP_ERROR_IO; } } else { + if (access(file_name, F_OK) != 0) { + return BSP_ERROR_IO; + } + bsp_prepare_hdf5_runtime(); hid_t f = H5Fopen(file_name, H5F_ACC_RDONLY, H5P_DEFAULT); + if (f < 0) { + return BSP_ERROR_IO; + } hid_t g = H5Gopen1(f, group); - bsp_matrix_t matrix = bsp_read_matrix_from_group(g); + if (g < 0) { + H5Fclose(f); + return BSP_ERROR_IO; + } + bsp_error_t error = + bsp_read_matrix_from_group_allocator(matrix, g, allocator); H5Gclose(g); H5Fclose(f); - return matrix; + return error; } } + +bsp_error_t bsp_read_matrix(bsp_matrix_t* matrix, const char* file_name, + const char* group) { + return bsp_read_matrix_allocator(matrix, file_name, group, + bsp_default_allocator); +} diff --git a/src/read_tensor.c b/src/read_tensor.c index 01453ac..cce6004 100644 --- a/src/read_tensor.c +++ b/src/read_tensor.c @@ -32,7 +32,11 @@ char* key_with_index(const char* key, size_t index) { bsp_tensor_t bsp_read_tensor_from_group(hid_t f) { bsp_tensor_t tensor = bsp_construct_default_tensor_t(); - char* json_string = bsp_read_attribute(f, (char*) "binsparse"); + char* json_string; + bsp_error_t error = bsp_read_attribute(&json_string, f, (char*) "binsparse"); + if (error != BSP_SUCCESS) { + return tensor; + } cJSON* j = cJSON_Parse(json_string); @@ -103,7 +107,12 @@ bsp_tensor_t bsp_read_tensor_from_group(hid_t f) { // base case: working with an element. if (strcmp(type, "element") == 0) { - bsp_array_t values = bsp_read_array(f, (char*) "values"); + bsp_array_t values; + bsp_error_t error = bsp_read_array(&values, f, (char*) "values"); + if (error != BSP_SUCCESS) { + free(json_string); + return tensor; + } cur_level->kind = BSP_TENSOR_ELEMENT; bsp_element_t* data = malloc(sizeof(bsp_element_t)); data->values = values; @@ -136,7 +145,12 @@ bsp_tensor_t bsp_read_tensor_from_group(hid_t f) { if (depth != 0) { char* pointers_key = key_with_index("pointers_to_", depth); data->pointers_to = malloc(sizeof(bsp_array_t)); - *data->pointers_to = bsp_read_array(f, pointers_key); + error = bsp_read_array(data->pointers_to, f, pointers_key); + if (error != BSP_SUCCESS) { + free(pointers_key); + free(json_string); + return tensor; + } free(pointers_key); } else { data->pointers_to = NULL; @@ -146,7 +160,12 @@ bsp_tensor_t bsp_read_tensor_from_group(hid_t f) { data->indices = malloc(rank * sizeof(bsp_array_t)); for (int idx = 0; idx < rank; idx++) { char* indices_key = key_with_index("indices_", depth + idx); - data->indices[idx] = bsp_read_array(f, indices_key); + error = bsp_read_array(&data->indices[idx], f, indices_key); + if (error != BSP_SUCCESS) { + free(indices_key); + free(json_string); + return tensor; + } free(indices_key); } diff --git a/src/write_matrix.c b/src/write_matrix.c index 7c012b0..38e38a6 100644 --- a/src/write_matrix.c +++ b/src/write_matrix.c @@ -5,11 +5,20 @@ */ #include -#include +#include +#include #include #include #include +static void bsp_prepare_hdf5_runtime(void) { + static bool initialized = false; + if (!initialized) { + H5dont_atexit(); + initialized = true; + } +} + char* bsp_generate_json(bsp_matrix_t matrix, cJSON* user_json) { cJSON* j = cJSON_CreateObject(); assert(j != NULL); @@ -20,9 +29,12 @@ char* bsp_generate_json(bsp_matrix_t matrix, cJSON* user_json) { cJSON_AddItemToObject(j, "binsparse", binsparse); - cJSON* item; - cJSON_ArrayForEach(item, user_json) { - cJSON_AddItemToObject(j, item->string, item); + if (user_json != NULL) { + cJSON* item; + cJSON_ArrayForEach(item, user_json) { + cJSON* item_copy = cJSON_Duplicate(item, 1); // 1 = deep copy + cJSON_AddItemToObject(j, item->string, item_copy); + } } cJSON_AddStringToObject(binsparse, "version", BINSPARSE_VERSION); @@ -86,63 +98,115 @@ char* bsp_generate_json(bsp_matrix_t matrix, cJSON* user_json) { return string; } -int bsp_write_matrix_to_group(hid_t f, bsp_matrix_t matrix, cJSON* user_json, - int compression_level) { - int result = +bsp_error_t bsp_write_matrix_to_group_cjson(hid_t f, bsp_matrix_t matrix, + cJSON* user_json, + int compression_level) { + bsp_error_t error = bsp_write_array(f, (char*) "values", matrix.values, compression_level); - - if (result != 0) - return result; + if (error != BSP_SUCCESS) { + return error; + } if (matrix.indices_0.size > 0) { - result = bsp_write_array(f, (char*) "indices_0", matrix.indices_0, - compression_level); - if (result != 0) { - return result; + error = bsp_write_array(f, (char*) "indices_0", matrix.indices_0, + compression_level); + if (error != BSP_SUCCESS) { + return error; } } if (matrix.indices_1.size > 0) { - result = bsp_write_array(f, (char*) "indices_1", matrix.indices_1, - compression_level); - if (result != 0) { - return result; + error = bsp_write_array(f, (char*) "indices_1", matrix.indices_1, + compression_level); + if (error != BSP_SUCCESS) { + return error; } } if (matrix.pointers_to_1.size > 0) { - result = bsp_write_array(f, (char*) "pointers_to_1", matrix.pointers_to_1, - compression_level); - if (result != 0) { - return result; + error = bsp_write_array(f, (char*) "pointers_to_1", matrix.pointers_to_1, + compression_level); + if (error != BSP_SUCCESS) { + return error; } } char* json_string = bsp_generate_json(matrix, user_json); - bsp_write_attribute(f, (char*) "binsparse", json_string); + error = bsp_write_attribute(f, (char*) "binsparse", json_string); + if (error != BSP_SUCCESS) { + free(json_string); + return error; + } free(json_string); - return 0; + return BSP_SUCCESS; +} + +bsp_error_t bsp_write_matrix_to_group(hid_t f, bsp_matrix_t matrix, + const char* user_json, + int compression_level) { + cJSON* user_json_cjson = NULL; + if (user_json != NULL) { + user_json_cjson = cJSON_Parse(user_json); + } + if (user_json_cjson == NULL) { + user_json_cjson = cJSON_CreateObject(); + } + bsp_error_t error = bsp_write_matrix_to_group_cjson( + f, matrix, user_json_cjson, compression_level); + cJSON_Delete(user_json_cjson); + return error; } -int bsp_write_matrix(const char* fname, bsp_matrix_t matrix, const char* group, - cJSON* user_json, int compression_level) { +bsp_error_t bsp_write_matrix_cjson(const char* fname, bsp_matrix_t matrix, + const char* group, cJSON* user_json, + int compression_level) { if (group == NULL) { + bsp_prepare_hdf5_runtime(); hid_t f = H5Fcreate(fname, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); - bsp_write_matrix_to_group(f, matrix, user_json, compression_level); + bsp_error_t error = bsp_write_matrix_to_group_cjson(f, matrix, user_json, + compression_level); + if (error != BSP_SUCCESS) { + H5Fclose(f); + return error; + } H5Fclose(f); } else { hid_t f; + bsp_prepare_hdf5_runtime(); if (access(fname, F_OK) == 0) { f = H5Fopen(fname, H5F_ACC_RDWR, H5P_DEFAULT); } else { f = H5Fcreate(fname, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); } hid_t g = H5Gcreate1(f, group, H5P_DEFAULT); - bsp_write_matrix_to_group(g, matrix, user_json, compression_level); + bsp_error_t error = bsp_write_matrix_to_group_cjson(g, matrix, user_json, + compression_level); + if (error != BSP_SUCCESS) { + H5Gclose(g); + H5Fclose(f); + return error; + } H5Gclose(g); H5Fclose(f); } - return 0; + return BSP_SUCCESS; +} + +bsp_error_t bsp_write_matrix(const char* fname, bsp_matrix_t matrix, + const char* group, const char* user_json, + int compression_level) { + cJSON* user_json_cjson = NULL; + if (user_json != NULL) { + user_json_cjson = cJSON_Parse(user_json); + } + if (user_json_cjson == NULL) { + user_json_cjson = cJSON_CreateObject(); + } + + bsp_error_t error = bsp_write_matrix_cjson( + fname, matrix, group, user_json_cjson, compression_level); + cJSON_Delete(user_json_cjson); + return error; } diff --git a/src/write_tensor.c b/src/write_tensor.c index 11a84ef..1b71f17 100644 --- a/src/write_tensor.c +++ b/src/write_tensor.c @@ -4,13 +4,15 @@ * SPDX-License-Identifier: BSD-3-Clause */ +#include +#include +#include + #include -#include #include -#include -#include #include +#include static cJSON* init_tensor_json(bsp_tensor_t tensor, cJSON* user_json) { cJSON* j = cJSON_CreateObject(); @@ -53,8 +55,8 @@ static cJSON* init_tensor_json(bsp_tensor_t tensor, cJSON* user_json) { return j; } -int bsp_write_tensor_to_group(hid_t f, bsp_tensor_t tensor, cJSON* user_json, - int compression_level) { +bsp_error_t bsp_write_tensor_to_group(hid_t f, bsp_tensor_t tensor, + cJSON* user_json, int compression_level) { // bsp_matrix_t matrix; cJSON* j = init_tensor_json(tensor, user_json); // tensor: @@ -85,10 +87,11 @@ int bsp_write_tensor_to_group(hid_t f, bsp_tensor_t tensor, cJSON* user_json, } // attempt to write an array. - int result = bsp_write_array(f, (char*) "values", values, compression_level); - if (result != 0) { + bsp_error_t error = + bsp_write_array(f, (char*) "values", values, compression_level); + if (error != BSP_SUCCESS) { cJSON_Delete(j); - return result; + return error; } int rank = 0; @@ -107,11 +110,11 @@ int bsp_write_tensor_to_group(hid_t f, bsp_tensor_t tensor, cJSON* user_json, cJSON_AddStringToObject(data_types, key_with_index("pointers_to_", rank), bsp_get_type_string(sparse->pointers_to->type)); - result = bsp_write_array(f, key_with_index("pointers_to_", rank), - *sparse->pointers_to, compression_level); - if (result != 0) { + error = bsp_write_array(f, key_with_index("pointers_to_", rank), + *sparse->pointers_to, compression_level); + if (error != BSP_SUCCESS) { cJSON_Delete(j); - return result; + return error; } } @@ -119,11 +122,11 @@ int bsp_write_tensor_to_group(hid_t f, bsp_tensor_t tensor, cJSON* user_json, cJSON_AddStringToObject(data_types, key_with_index("indices_", rank + i), bsp_get_type_string(sparse->indices[i].type)); - result = bsp_write_array(f, key_with_index("indices_", rank + i), - sparse->indices[i], compression_level); - if (result != 0) { + error = bsp_write_array(f, key_with_index("indices_", rank + i), + sparse->indices[i], compression_level); + if (error != BSP_SUCCESS) { cJSON_Delete(j); - return result; + return error; } } @@ -153,17 +156,28 @@ int bsp_write_tensor_to_group(hid_t f, bsp_tensor_t tensor, cJSON* user_json, } char* json_string = cJSON_Print(j); - bsp_write_attribute(f, (char*) "binsparse", json_string); + error = bsp_write_attribute(f, (char*) "binsparse", json_string); + if (error != BSP_SUCCESS) { + free(json_string); + cJSON_Delete(j); + return error; + } free(json_string); - return 0; + return BSP_SUCCESS; } -int bsp_write_tensor(const char* fname, bsp_tensor_t tensor, const char* group, - cJSON* user_json, int compression_level) { +bsp_error_t bsp_write_tensor(const char* fname, bsp_tensor_t tensor, + const char* group, cJSON* user_json, + int compression_level) { if (group == NULL) { hid_t f = H5Fcreate(fname, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); - bsp_write_tensor_to_group(f, tensor, user_json, compression_level); + bsp_error_t error = + bsp_write_tensor_to_group(f, tensor, user_json, compression_level); + if (error != BSP_SUCCESS) { + H5Fclose(f); + return error; + } H5Fclose(f); } else { hid_t f; @@ -177,5 +191,5 @@ int bsp_write_tensor(const char* fname, bsp_tensor_t tensor, const char* group, H5Gclose(g); H5Fclose(f); } - return 0; + return BSP_SUCCESS; }