diff --git a/bindings/matlab/Contents.m b/bindings/matlab/Contents.m index c785ccc..27ef4ce 100644 --- a/bindings/matlab/Contents.m +++ b/bindings/matlab/Contents.m @@ -6,6 +6,7 @@ % 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 +% binsparse_write_ssmc_coo - write MATLAB CSC directly as row-sorted COO % binsparse_read_string_dataset - read an HDF5 UTF-8 string dataset % % MATLAB helpers: @@ -24,9 +25,11 @@ % 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_binsparse_matrix_views - non-owning matrix view ownership tests % test_binsparse_struct - tests for the Binsparse struct helpers % test_binsparse_to_ssmc_problem - tests for Binsparse Problem conversion % test_binsparse_write_ssmc_problem - end-to-end SSMC writer test +% test_binsparse_write_ssmc_coo - direct row-sorted COO writer tests % test_binsparse_roundtrip_dir - round-trip every .h5 file in a directory % SPDX-FileCopyrightText: 2024 Binsparse Developers diff --git a/bindings/matlab/README.md b/bindings/matlab/README.md index b0aee5d..e3292a3 100644 --- a/bindings/matlab/README.md +++ b/bindings/matlab/README.md @@ -186,6 +186,7 @@ end | `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_ssmc_coo.c` | Direct memory-efficient MATLAB CSC to row-sorted COO writer | | `binsparse_write_string_dataset.c` | MEX function writing HDF5 UTF-8 string datasets | | `matlab_bsp_helpers.h` | Shared MATLAB/Binsparse conversion helpers for the MEX sources | | `binsparse_write_ssmc_problem.m` | Write a full SSMC Problem struct to one Binsparse file | @@ -199,9 +200,11 @@ end | `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_binsparse_matrix_views.m` | Test non-owning matrix view ownership and cleanup | | `test_binsparse_struct.m` | Test script for the struct helpers | | `test_binsparse_to_ssmc_problem.m` | Test script for Problem conversion | | `test_binsparse_write_ssmc_problem.m` | End-to-end test for the SSMC writer | +| `test_binsparse_write_ssmc_coo.m` | Test direct row-sorted COO conversion and writing | | `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 | diff --git a/bindings/matlab/binsparse_build_matlab_bindings.m b/bindings/matlab/binsparse_build_matlab_bindings.m index 701d837..4347498 100644 --- a/bindings/matlab/binsparse_build_matlab_bindings.m +++ b/bindings/matlab/binsparse_build_matlab_bindings.m @@ -107,6 +107,7 @@ function compile_mex_functions(paths, verbose) % List of MEX functions to compile mex_files = {'binsparse_read.c', 'binsparse_write.c', ... 'binsparse_from_ssmc.c', 'binsparse_minimize_types.c', ... + 'binsparse_write_ssmc_coo.c', ... 'binsparse_write_string_dataset.c', ... 'binsparse_read_string_dataset.c'}; @@ -140,7 +141,8 @@ function compile_mex_functions(paths, verbose) rpath = [' LDFLAGS=''$LDFLAGS -fPIC ' rpath ' '' '] ; end - include_flags = sprintf('-I%s', paths.include_dir); + include_flags = sprintf('-I%s -I%s', paths.include_dir, ... + fullfile(lib_dir, 'include')); if ~isempty(paths.hdf5_include_dir) include_flags = sprintf('%s -I%s', include_flags, ... paths.hdf5_include_dir); diff --git a/bindings/matlab/binsparse_build_octave_bindings.m b/bindings/matlab/binsparse_build_octave_bindings.m index c9ba44d..cd85aee 100644 --- a/bindings/matlab/binsparse_build_octave_bindings.m +++ b/bindings/matlab/binsparse_build_octave_bindings.m @@ -120,6 +120,7 @@ function compile_octave_functions(paths, verbose) % List of MEX functions to compile mex_files = {'binsparse_read.c', 'binsparse_write.c', ... 'binsparse_from_ssmc.c', 'binsparse_minimize_types.c', ... + 'binsparse_write_ssmc_coo.c', ... 'binsparse_write_string_dataset.c', ... 'binsparse_read_string_dataset.c'}; @@ -136,12 +137,13 @@ function compile_octave_functions(paths, verbose) fprintf(' Compiling %s... ', mex_file); % Prepare mkoctfile command with library linking - include_flag = sprintf('-I%s', paths.include_dir); + lib_dir = fullfile(paths.binsparse_root, 'build'); + include_flag = sprintf('-I%s -I%s', paths.include_dir, ... + fullfile(lib_dir, 'include')); 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'); diff --git a/bindings/matlab/binsparse_from_ssmc.c b/bindings/matlab/binsparse_from_ssmc.c index c91722e..03259f8 100644 --- a/bindings/matlab/binsparse_from_ssmc.c +++ b/bindings/matlab/binsparse_from_ssmc.c @@ -64,6 +64,46 @@ static bsp_type_t sparse_value_type(const matlab_csc_t* matrix) { return matrix->is_complex ? BSP_COMPLEX_FLOAT64 : BSP_FLOAT64; } +static bsp_type_t index_type_for_extent(size_t extent) { + return bsp_pick_integer_type(extent == 0 ? 0 : extent - 1); +} + +static bool sparse_values_are_iso(const matlab_csc_t* matrix, double* real, + double* imag) { + if (matrix->nnz == 0) { + return false; + } + + *real = matrix->values[0]; + *imag = matrix->imag_values ? matrix->imag_values[0] : 0.0; + for (size_t i = 1; i < matrix->nnz; i++) { + double current_imag = matrix->imag_values ? matrix->imag_values[i] : 0.0; + if (matrix->values[i] != *real || current_imag != *imag) { + return false; + } + } + return true; +} + +static bool merged_values_are_iso(const matlab_csc_t* a, + const matlab_csc_t* z, double* real, + double* imag) { + if (a->nnz + z->nnz == 0) { + return false; + } + + if (a->nnz == 0) { + *real = 0.0; + *imag = 0.0; + return true; + } + + if (!sparse_values_are_iso(a, real, imag)) { + return false; + } + return z->nnz == 0 || (*real == 0.0 && *imag == 0.0); +} + static void write_sparse_value(const matlab_csc_t* matrix, mwIndex src, bsp_matrix_t* out, uint64_t dst) { if (matrix->is_complex) { @@ -86,6 +126,16 @@ static void write_explicit_zero(bsp_matrix_t* out, uint64_t dst) { } } +static void write_iso_value(bsp_matrix_t* out, double real, double imag) { + if (out->values.type == BSP_COMPLEX_FLOAT64) { + double _Complex* out_values = (double _Complex*) out->values.data; + out_values[0] = real + imag * I; + } else { + double* out_values = (double*) out->values.data; + out_values[0] = real; + } +} + static void build_csc_merged(const matlab_csc_t* a, const matlab_csc_t* z, bsp_matrix_t* out) { bsp_error_t error; @@ -96,37 +146,45 @@ static void build_csc_merged(const matlab_csc_t* a, const matlab_csc_t* z, out->nnz = a->nnz + z->nnz; out->format = BSP_CSC; out->structure = BSP_GENERAL; - out->is_iso = false; + double iso_real = 0.0; + double iso_imag = 0.0; + out->is_iso = merged_values_are_iso(a, z, &iso_real, &iso_imag); error = construct_array_with_allocator( - &out->values, out->nnz, sparse_value_type(a), bsp_matlab_allocator); + &out->values, out->is_iso ? 1 : 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); + bsp_type_t row_index_type = index_type_for_extent(out->nrows); + error = construct_array_with_allocator(&out->indices_1, out->nnz, + row_index_type, bsp_matlab_allocator); if (error != BSP_SUCCESS) { mexErrMsgIdAndTxt("BinSparse:MemoryError", "Failed to allocate indices array"); } + bsp_type_t pointer_type = bsp_pick_integer_type(out->nnz); error = construct_array_with_allocator(&out->pointers_to_1, out->ncols + 1, - BSP_UINT64, bsp_matlab_allocator); + pointer_type, 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; + if (out->is_iso) { + write_iso_value(out, iso_real, iso_imag); + } - out_colptr[0] = 0; + uint64_t running = 0; + bsp_array_write(out->pointers_to_1, 0, running); 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; + running += a_count + z_count; + bsp_array_write(out->pointers_to_1, j + 1, running); } for (mwIndex j = 0; j < a->ncols; j++) { @@ -134,18 +192,22 @@ static void build_csc_merged(const matlab_csc_t* a, const matlab_csc_t* z, 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]; + uint64_t out_ptr = (uint64_t) a->colptr[j] + (uint64_t) z->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]; + if (!out->is_iso) { + write_sparse_value(a, a_ptr, out, out_ptr); + } + bsp_array_write(out->indices_1, out_ptr, 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]; + if (!out->is_iso) { + write_explicit_zero(out, out_ptr); + } + bsp_array_write(out->indices_1, out_ptr, z->rowind[z_ptr]); z_ptr++; } else { mexErrMsgIdAndTxt("BinSparse:DuplicateIndex", @@ -154,7 +216,8 @@ static void build_csc_merged(const matlab_csc_t* a, const matlab_csc_t* z, out_ptr++; } - if (out_ptr != out_colptr[j + 1]) { + if (out_ptr != (uint64_t) a->colptr[j + 1] + + (uint64_t) z->colptr[j + 1]) { mexErrMsgIdAndTxt("BinSparse:InternalError", "Merged column counts do not match"); } @@ -204,39 +267,47 @@ static void build_csc_from_a(const matlab_csc_t* a, bsp_matrix_t* out) { out->nnz = a->nnz; out->format = BSP_CSC; out->structure = BSP_GENERAL; - out->is_iso = false; + double iso_real = 0.0; + double iso_imag = 0.0; + out->is_iso = sparse_values_are_iso(a, &iso_real, &iso_imag); error = construct_array_with_allocator( - &out->values, out->nnz, sparse_value_type(a), bsp_matlab_allocator); + &out->values, out->is_iso ? 1 : 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); + bsp_type_t row_index_type = index_type_for_extent(out->nrows); + error = construct_array_with_allocator(&out->indices_1, out->nnz, + row_index_type, bsp_matlab_allocator); if (error != BSP_SUCCESS) { mexErrMsgIdAndTxt("BinSparse:MemoryError", "Failed to allocate indices array"); } + bsp_type_t pointer_type = bsp_pick_integer_type(out->nnz); error = construct_array_with_allocator(&out->pointers_to_1, out->ncols + 1, - BSP_UINT64, bsp_matlab_allocator); + pointer_type, 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; + if (out->is_iso) { + write_iso_value(out, iso_real, iso_imag); + } 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]; + if (!out->is_iso) { + write_sparse_value(a, i, out, (uint64_t) i); + } + bsp_array_write(out->indices_1, i, a->rowind[i]); } for (size_t i = 0; i < out->ncols + 1; i++) { - out_colptr[i] = (uint64_t) a->colptr[i]; + bsp_array_write(out->pointers_to_1, i, a->colptr[i]); } } diff --git a/bindings/matlab/binsparse_minimize_types.c b/bindings/matlab/binsparse_minimize_types.c index 04ef621..df075e3 100644 --- a/bindings/matlab/binsparse_minimize_types.c +++ b/bindings/matlab/binsparse_minimize_types.c @@ -429,8 +429,7 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { } bsp_matrix_t matrix; - bsp_error_t error = matlab_struct_to_bsp_matrix_allocator( - prhs[0], &matrix, bsp_matlab_allocator); + bsp_error_t error = matlab_struct_to_bsp_matrix_view(prhs[0], &matrix); if (error != BSP_SUCCESS) { mexErrMsgIdAndTxt("BinSparse:ConversionError", "Failed to convert MATLAB struct to matrix: %s", diff --git a/bindings/matlab/binsparse_write.c b/bindings/matlab/binsparse_write.c index 53e39b3..3ec0e6e 100644 --- a/bindings/matlab/binsparse_write.c +++ b/bindings/matlab/binsparse_write.c @@ -78,8 +78,7 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { } // Convert MATLAB struct to bsp_matrix_t - error = matlab_struct_to_bsp_matrix_allocator(prhs[1], &matrix, - bsp_matlab_allocator); + error = matlab_struct_to_bsp_matrix_view(prhs[1], &matrix); if (error != BSP_SUCCESS) { mxFree(filename); mexErrMsgIdAndTxt("BinSparse:ConversionError", diff --git a/bindings/matlab/binsparse_write_ssmc_coo.c b/bindings/matlab/binsparse_write_ssmc_coo.c new file mode 100644 index 0000000..542716d --- /dev/null +++ b/bindings/matlab/binsparse_write_ssmc_coo.c @@ -0,0 +1,504 @@ +/* + * SPDX-FileCopyrightText: 2026 Binsparse Developers + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +/** + * Write a MATLAB CSC sparse matrix directly as row-sorted Binsparse COO. + * + * Usage: + * binsparse_write_ssmc_coo(filename, A, Zeros, group, json, compression) + * + * The conversion uses a counting scatter rather than the generic COO qsort. + * Only reordered columns and non-ISO values are materialized. Row indices are + * generated directly into bounded HDF5 write buffers from the row endpoints. + */ + +#include "matlab_bsp_helpers.h" +#include "matlab_bsp_strings.h" +#include "mex.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#define COO_WRITE_CHUNK_BYTES ((size_t) (1024 * 1024)) + +typedef struct value_analysis_t { + bool is_iso; + bool fits_float32; + double iso_real; + double iso_imag; +} value_analysis_t; + +static bsp_type_t index_type_for_extent(size_t extent) { + return bsp_pick_integer_type(extent == 0 ? 0 : extent - 1); +} + +static char* required_string(const mxArray* value, const char* name, + bool allow_empty) { + if (!mxIsChar(value) || mxGetM(value) > 1 || + (!allow_empty && mxIsEmpty(value))) { + mexErrMsgIdAndTxt("BinSparse:InvalidString", + "%s must be a character row vector", name); + } + char* string = mxArrayToString(value); + if (!string) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", "Failed to read %s", name); + } + return string; +} + +static int compression_level(const mxArray* value) { + if (!mxIsNumeric(value) || mxIsComplex(value) || + mxGetNumberOfElements(value) != 1) { + mexErrMsgIdAndTxt("BinSparse:InvalidCompression", + "Compression level must be a real numeric scalar"); + } + double level = mxGetScalar(value); + if (level < 0.0 || level > 9.0 || level != (double) (int) level) { + mexErrMsgIdAndTxt("BinSparse:InvalidCompression", + "Compression level must be an integer from 0 to 9"); + } + return (int) level; +} + +static bool value_fits_float32(double real, double imag, bool is_complex) { + if ((double) (float) real != real) { + return false; + } + return !is_complex || (double) (float) imag == imag; +} + +static value_analysis_t analyze_values(const matlab_csc_t* a, + const matlab_csc_t* z) { + value_analysis_t analysis = {.is_iso = a->nnz + z->nnz > 0, + .fits_float32 = true, + .iso_real = a->nnz > 0 ? a->values[0] : 0.0, + .iso_imag = + a->nnz > 0 && a->imag_values + ? a->imag_values[0] + : 0.0}; + + for (size_t i = 0; i < a->nnz; i++) { + double real = a->values[i]; + double imag = a->imag_values ? a->imag_values[i] : 0.0; + if (i > 0 && + (real != analysis.iso_real || imag != analysis.iso_imag)) { + analysis.is_iso = false; + } + if (!value_fits_float32(real, imag, a->is_complex)) { + analysis.fits_float32 = false; + } + } + + if (z->nnz > 0 && + (analysis.iso_real != 0.0 || analysis.iso_imag != 0.0)) { + analysis.is_iso = false; + } + return analysis; +} + +static bsp_type_t output_value_type(const matlab_csc_t* a, + value_analysis_t analysis) { + if (a->is_complex) { + return analysis.fits_float32 ? BSP_COMPLEX_FLOAT32 + : BSP_COMPLEX_FLOAT64; + } + return analysis.fits_float32 ? BSP_FLOAT32 : BSP_FLOAT64; +} + +static void write_output_value(bsp_array_t values, size_t index, double real, + double imag) { + switch (values.type) { + case BSP_FLOAT32: + ((float*) values.data)[index] = (float) real; + break; + case BSP_FLOAT64: + ((double*) values.data)[index] = real; + break; + case BSP_COMPLEX_FLOAT32: + ((float _Complex*) values.data)[index] = + (float) real + (float) imag * I; + break; + case BSP_COMPLEX_FLOAT64: + ((double _Complex*) values.data)[index] = real + imag * I; + break; + default: + break; + } +} + +static void write_column(bsp_array_t columns, size_t index, size_t column) { + bsp_array_write(columns, index, column); +} + +static void fill_indices(void* data, bsp_type_t type, size_t offset, + size_t count, size_t value) { + switch (type) { + case BSP_UINT8: { + uint8_t* x = (uint8_t*) data; + for (size_t k = 0; k < count; k++) + x[offset + k] = (uint8_t) value; + break; + } + case BSP_UINT16: { + uint16_t* x = (uint16_t*) data; + for (size_t k = 0; k < count; k++) + x[offset + k] = (uint16_t) value; + break; + } + case BSP_UINT32: { + uint32_t* x = (uint32_t*) data; + for (size_t k = 0; k < count; k++) + x[offset + k] = (uint32_t) value; + break; + } + case BSP_UINT64: { + uint64_t* x = (uint64_t*) data; + for (size_t k = 0; k < count; k++) + x[offset + k] = (uint64_t) value; + break; + } + default: + break; + } +} + +static bsp_error_t write_generated_rows(hid_t parent, const uint64_t* row_ends, + size_t nrows, size_t nnz, + bsp_type_t index_type, + int compression) { + hsize_t dimensions[1] = {(hsize_t) nnz}; + hid_t file_space = H5Screate_simple(1, dimensions, NULL); + if (file_space == H5I_INVALID_HID) { + return BSP_ERROR_IO; + } + + hid_t properties = H5Pcreate(H5P_DATASET_CREATE); + if (properties == H5I_INVALID_HID) { + H5Sclose(file_space); + return BSP_ERROR_IO; + } + + size_t element_size = bsp_type_size(index_type); + size_t capacity = COO_WRITE_CHUNK_BYTES / element_size; + if (capacity == 0) + capacity = 1; + if (nnz > 0 && capacity > nnz) + capacity = nnz; + + if (nnz > 0) { + hsize_t chunk[1] = {(hsize_t) capacity}; + if (H5Pset_chunk(properties, 1, chunk) < 0 || + (compression > 0 && + H5Pset_deflate(properties, (unsigned) compression) < 0)) { + H5Pclose(properties); + H5Sclose(file_space); + return BSP_ERROR_IO; + } + } + + hid_t dataset = H5Dcreate2(parent, "indices_0", + bsp_get_hdf5_standard_type(index_type), file_space, + H5P_DEFAULT, properties, H5P_DEFAULT); + H5Pclose(properties); + if (dataset == H5I_INVALID_HID) { + H5Sclose(file_space); + return BSP_ERROR_IO; + } + + if (nnz == 0) { + H5Dclose(dataset); + H5Sclose(file_space); + return BSP_SUCCESS; + } + + void* buffer = mxMalloc(capacity * element_size); + if (!buffer) { + H5Dclose(dataset); + H5Sclose(file_space); + return BSP_ERROR_MEMORY; + } + + size_t row = 0; + size_t written = 0; + bsp_error_t error = BSP_SUCCESS; + while (written < nnz) { + size_t buffered = 0; + while (buffered < capacity && written + buffered < nnz) { + while (row < nrows && written + buffered >= row_ends[row]) + row++; + if (row == nrows) { + error = BSP_ERROR_INTERNAL; + break; + } + size_t available = (size_t) row_ends[row] - (written + buffered); + size_t room = capacity - buffered; + size_t count = available < room ? available : room; + fill_indices(buffer, index_type, buffered, count, row); + buffered += count; + } + if (error != BSP_SUCCESS) + break; + + hsize_t start[1] = {(hsize_t) written}; + hsize_t count[1] = {(hsize_t) buffered}; + hid_t memory_space = H5Screate_simple(1, count, NULL); + if (memory_space == H5I_INVALID_HID || + H5Sselect_hyperslab(file_space, H5S_SELECT_SET, start, NULL, count, + NULL) < 0 || + H5Dwrite(dataset, bsp_get_hdf5_native_type(index_type), memory_space, + file_space, H5P_DEFAULT, buffer) < 0) { + if (memory_space != H5I_INVALID_HID) + H5Sclose(memory_space); + error = BSP_ERROR_IO; + break; + } + H5Sclose(memory_space); + written += buffered; + } + + mxFree(buffer); + H5Dclose(dataset); + H5Sclose(file_space); + return error; +} + +static bsp_error_t write_descriptor(hid_t parent, bsp_matrix_t descriptor, + const char* json) { + cJSON* user_json = json ? cJSON_Parse(json) : NULL; + if (!user_json) + user_json = cJSON_CreateObject(); + if (!user_json) + return BSP_ERROR_MEMORY; + + char* descriptor_json = bsp_generate_json(descriptor, user_json); + cJSON_Delete(user_json); + if (!descriptor_json) + return BSP_ERROR_MEMORY; + bsp_error_t error = + bsp_write_attribute(parent, (char*) "binsparse", descriptor_json); + free(descriptor_json); + return error; +} + +static bsp_error_t write_coo_file(const char* filename, const char* group, + bsp_matrix_t matrix, + const uint64_t* row_ends, const char* json, + int compression) { + hid_t file = H5I_INVALID_HID; + hid_t parent = H5I_INVALID_HID; + bool have_group = group && group[0] != '\0'; + + H5dont_atexit(); + if (!have_group) { + file = H5Fcreate(filename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + parent = file; + } else { + file = access(filename, F_OK) == 0 + ? H5Fopen(filename, H5F_ACC_RDWR, H5P_DEFAULT) + : H5Fcreate(filename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + if (file != H5I_INVALID_HID) + parent = H5Gcreate1(file, group, H5P_DEFAULT); + } + if (file == H5I_INVALID_HID || parent == H5I_INVALID_HID) { + if (file != H5I_INVALID_HID) + H5Fclose(file); + return BSP_ERROR_IO; + } + + bsp_error_t error = + bsp_write_array(parent, "values", matrix.values, compression); + if (error == BSP_SUCCESS) { + error = write_generated_rows(parent, row_ends, matrix.nrows, matrix.nnz, + matrix.indices_0.type, compression); + } + if (error == BSP_SUCCESS) { + error = bsp_write_array(parent, "indices_1", matrix.indices_1, compression); + } + if (error == BSP_SUCCESS) { + error = write_descriptor(parent, matrix, json); + } + + if (have_group) + H5Gclose(parent); + H5Fclose(file); + return error; +} + +void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { + (void) plhs; + bsp_lock_mex_module(); + if (nrhs != 6) { + mexErrMsgIdAndTxt( + "BinSparse:InvalidArgs", + "Usage: binsparse_write_ssmc_coo(filename, A, Zeros, group, json, " + "compression)"); + } + if (nlhs != 0) { + mexErrMsgIdAndTxt("BinSparse:TooManyOutputs", "No output is returned"); + } + + char* filename = required_string(prhs[0], "filename", false); + char* group = required_string(prhs[3], "group", true); + char* json = NULL; + if (!mxIsEmpty(prhs[4])) { + json = bsp_mx_to_utf8(prhs[4]); + if (!json) { + mxFree(group); + mxFree(filename); + mexErrMsgIdAndTxt("BinSparse:InvalidJSON", + "JSON must be a valid character vector"); + } + } + int compression = compression_level(prhs[5]); + + const mxArray* mx_a = prhs[1]; + const mxArray* mx_z = mxIsEmpty(prhs[2]) ? NULL : prhs[2]; + if (!mxIsSparse(mx_a) || !mxIsDouble(mx_a)) { + mexErrMsgIdAndTxt("BinSparse:InvalidMatrix", + "A must be a sparse double matrix"); + } + if (mx_z && + (!mxIsSparse(mx_z) || mxIsComplex(mx_z) || + (!mxIsDouble(mx_z) && !mxIsLogical(mx_z)))) { + mexErrMsgIdAndTxt("BinSparse:InvalidZeros", + "Zeros must be a real sparse double or logical matrix"); + } + if (mx_z && + (mxGetM(mx_a) != mxGetM(mx_z) || mxGetN(mx_a) != mxGetN(mx_z))) { + mexErrMsgIdAndTxt("BinSparse:DimensionMismatch", + "A and Zeros must have matching dimensions"); + } + + matlab_csc_t a = {0}; + matlab_csc_t z = {0}; + if (extract_matlab_csc(mx_a, &a) != 0 || + (mx_z && extract_matlab_csc(mx_z, &z) != 0)) { + mexErrMsgIdAndTxt("BinSparse:InvalidMatrix", + "Failed to access MATLAB CSC arrays"); + } + if (a.nnz > SIZE_MAX - z.nnz) { + mexErrMsgIdAndTxt("BinSparse:SizeOverflow", "Stored entry count overflows"); + } + size_t nnz = a.nnz + z.nnz; + + uint64_t* positions = + a.nrows > 0 ? (uint64_t*) mxCalloc(a.nrows, sizeof(uint64_t)) : NULL; + if (a.nrows > 0 && !positions) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to allocate row counts"); + } + for (size_t p = 0; p < a.nnz; p++) + positions[a.rowind[p]]++; + for (size_t p = 0; p < z.nnz; p++) + positions[z.rowind[p]]++; + + uint64_t running = 0; + for (size_t row = 0; row < a.nrows; row++) { + uint64_t count = positions[row]; + positions[row] = running; + running += count; + } + if (running != nnz) { + mxFree(positions); + mexErrMsgIdAndTxt("BinSparse:InternalError", "Row counts do not match nnz"); + } + + value_analysis_t analysis = analyze_values(&a, &z); + bsp_matrix_t output; + bsp_construct_default_matrix_t_allocator(&output, bsp_matlab_allocator); + output.nrows = a.nrows; + output.ncols = a.ncols; + output.nnz = nnz; + output.is_iso = analysis.is_iso; + output.format = BSP_COOR; + output.structure = BSP_GENERAL; + + bsp_type_t value_type = output_value_type(&a, analysis); + bsp_error_t error = bsp_construct_array_t_allocator( + &output.values, output.is_iso ? 1 : nnz, value_type, + bsp_matlab_allocator); + if (error == BSP_SUCCESS) { + error = bsp_construct_array_t_allocator( + &output.indices_1, nnz, index_type_for_extent(a.ncols), + bsp_matlab_allocator); + } + if (error != BSP_SUCCESS) { + bsp_destroy_matrix_t(&output); + mxFree(positions); + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to allocate COO output arrays"); + } + output.indices_0.size = nnz; + output.indices_0.type = index_type_for_extent(a.nrows); + output.indices_0.data = nnz > 0 ? output.indices_1.data : NULL; + output.indices_0.allocator = bsp_matlab_view_allocator; + + if (output.is_iso) { + write_output_value(output.values, 0, analysis.iso_real, analysis.iso_imag); + } + + const char* conversion_error = NULL; + for (size_t column = 0; column < a.ncols && !conversion_error; column++) { + mwIndex ap = a.colptr[column]; + mwIndex aend = a.colptr[column + 1]; + mwIndex zp = z.colptr ? z.colptr[column] : 0; + mwIndex zend = z.colptr ? z.colptr[column + 1] : 0; + while (ap < aend || zp < zend) { + bool take_a = zp >= zend || + (ap < aend && a.rowind[ap] < z.rowind[zp]); + bool take_z = ap >= aend || + (zp < zend && z.rowind[zp] < a.rowind[ap]); + if (!take_a && !take_z) { + conversion_error = "Duplicate indices between A and Zeros"; + break; + } + + size_t row = take_a ? a.rowind[ap] : z.rowind[zp]; + size_t destination = (size_t) positions[row]++; + write_column(output.indices_1, destination, column); + if (!output.is_iso) { + double real = take_a ? a.values[ap] : 0.0; + double imag = take_a && a.imag_values ? a.imag_values[ap] : 0.0; + write_output_value(output.values, destination, real, imag); + } + if (take_a) + ap++; + else + zp++; + } + } + + if (conversion_error) { + bsp_destroy_array_t(&output.values); + bsp_destroy_array_t(&output.indices_1); + mxFree(positions); + mxFree(json); + mxFree(group); + mxFree(filename); + mexErrMsgIdAndTxt("BinSparse:DuplicateIndex", "%s", conversion_error); + } + + error = write_coo_file(filename, group, output, positions, json, compression); + bsp_destroy_array_t(&output.values); + bsp_destroy_array_t(&output.indices_1); + mxFree(positions); + if (json) + mxFree(json); + mxFree(group); + mxFree(filename); + + if (error != BSP_SUCCESS) { + mexErrMsgIdAndTxt("BinSparse:WriteError", "Failed to write COO matrix: %s", + bsp_get_error_string(error)); + } +} diff --git a/bindings/matlab/binsparse_write_ssmc_problem.m b/bindings/matlab/binsparse_write_ssmc_problem.m index 41e8f4e..fa158dd 100644 --- a/bindings/matlab/binsparse_write_ssmc_problem.m +++ b/bindings/matlab/binsparse_write_ssmc_problem.m @@ -73,7 +73,20 @@ function binsparse_write_ssmc_problem(problem, output_filename, format, compress % Primary matrix A = P.A; -if issparse(A) +primary_metadata = metadata_json(P, 'A'); +if issparse(A) && any(strcmp(format, {'COO', 'COOR'})) + if exist('binsparse_write_ssmc_coo', 'file') ~= 3 + error('binsparse_write_ssmc_problem:MissingCOOWriter', ... + 'COO output requires binsparse_write_ssmc_coo on the path'); + end + if isfield(P, 'Zeros') && ~isempty(P.Zeros) + Zeros = P.Zeros; + else + Zeros = []; + end + binsparse_write_ssmc_coo(output_filename, A, Zeros, '', ... + primary_metadata, compression_level); +elseif issparse(A) if isfield(P, 'Zeros') && ~isempty(P.Zeros) Zeros = P.Zeros; mat = binsparse_from_ssmc(A, Zeros, format); @@ -84,9 +97,10 @@ function binsparse_write_ssmc_problem(problem, output_filename, format, compress 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); +if ~(issparse(A) && any(strcmp(format, {'COO', 'COOR'}))) + mat = binsparse_minimize_types(mat); + binsparse_write(output_filename, mat, '', primary_metadata, compression_level); +end % Handle aux struct if isfield(P, 'aux') && isstruct(P.aux) @@ -132,7 +146,11 @@ function handle_aux_entry(name, value, output_filename, format, compression_leve value = double(value); end - if issparse(value) + if issparse(value) && any(strcmp(format, {'COO', 'COOR'})) + binsparse_write_ssmc_coo(output_filename, value, [], name, ... + entry_metadata_json(name), compression_level); + return; + elseif issparse(value) bsp = binsparse_from_ssmc(value, format); elseif isnumeric(value) bsp = binsparse_from_ssmc(value, dense_format_for(value)); diff --git a/bindings/matlab/compile_octave.sh b/bindings/matlab/compile_octave.sh index 53e7f80..2826f6e 100755 --- a/bindings/matlab/compile_octave.sh +++ b/bindings/matlab/compile_octave.sh @@ -148,7 +148,7 @@ if [ "$CLEAN" = true ]; then fi # List of MEX files to compile -MEX_FILES=("binsparse_read.c" "binsparse_write.c" "binsparse_from_ssmc.c" "binsparse_minimize_types.c" "binsparse_write_string_dataset.c" "binsparse_read_string_dataset.c") +MEX_FILES=("binsparse_read.c" "binsparse_write.c" "binsparse_from_ssmc.c" "binsparse_minimize_types.c" "binsparse_write_ssmc_coo.c" "binsparse_write_string_dataset.c" "binsparse_read_string_dataset.c") print_info "Compiling MEX functions..." @@ -166,7 +166,7 @@ for mex_file in "${MEX_FILES[@]}"; do LIB_PATH="$LIB_DIR/libbinsparse.a" CJSON_LIB_DIR="$LIB_DIR/_deps/cjson-build" - INCLUDE_FLAGS="-I\"$INCLUDE_DIR\"" + INCLUDE_FLAGS="-I\"$INCLUDE_DIR\" -I\"$LIB_DIR/include\"" if [ -n "$HDF5_INCLUDE_DIR" ]; then INCLUDE_FLAGS="$INCLUDE_FLAGS -I\"$HDF5_INCLUDE_DIR\"" fi diff --git a/bindings/matlab/matlab_bsp_helpers.h b/bindings/matlab/matlab_bsp_helpers.h index 02c607f..1631a20 100644 --- a/bindings/matlab/matlab_bsp_helpers.h +++ b/bindings/matlab/matlab_bsp_helpers.h @@ -26,6 +26,20 @@ typedef struct { static const bsp_allocator_t bsp_matlab_allocator = {.malloc = mxMalloc, .free = mxFree}; +// Arrays using this allocator borrow their data from a MATLAB input array. +// They must never allocate or free through the Binsparse array API. Returning +// NULL from malloc makes accidental attempts to allocate through a view fail +// instead of creating memory that the no-op free would leak. +static inline void* bsp_matlab_view_malloc(size_t size) { + (void) size; + return NULL; +} + +static inline void bsp_matlab_view_free(void* data) { (void) data; } + +static const bsp_allocator_t bsp_matlab_view_allocator = { + .malloc = bsp_matlab_view_malloc, .free = bsp_matlab_view_free}; + static inline int extract_matlab_csc(const mxArray* mx_matrix, matlab_csc_t* csc_matrix) { if (!mx_matrix || !csc_matrix) { @@ -71,76 +85,91 @@ static inline int extract_matlab_csc(const mxArray* mx_matrix, 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); +static inline bsp_error_t matlab_bsp_array_type(const mxArray* mx_array, + bsp_type_t* bsp_type, + size_t* element_size) { + if (!mx_array || !bsp_type || !element_size || mxIsSparse(mx_array)) { + return BSP_ERROR_TYPE; + } - 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); + *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); + *bsp_type = BSP_COMPLEX_FLOAT32; + *element_size = sizeof(float _Complex); } else { - return BSP_INVALID_TYPE; + return BSP_ERROR_TYPE; } } else { switch (class_id) { case mxDOUBLE_CLASS: - bsp_type = BSP_FLOAT64; - element_size = sizeof(double); + *bsp_type = BSP_FLOAT64; + *element_size = sizeof(double); break; case mxSINGLE_CLASS: - bsp_type = BSP_FLOAT32; - element_size = sizeof(float); + *bsp_type = BSP_FLOAT32; + *element_size = sizeof(float); break; case mxUINT64_CLASS: - bsp_type = BSP_UINT64; - element_size = sizeof(uint64_t); + *bsp_type = BSP_UINT64; + *element_size = sizeof(uint64_t); break; case mxUINT32_CLASS: - bsp_type = BSP_UINT32; - element_size = sizeof(uint32_t); + *bsp_type = BSP_UINT32; + *element_size = sizeof(uint32_t); break; case mxUINT16_CLASS: - bsp_type = BSP_UINT16; - element_size = sizeof(uint16_t); + *bsp_type = BSP_UINT16; + *element_size = sizeof(uint16_t); break; case mxUINT8_CLASS: - bsp_type = BSP_UINT8; - element_size = sizeof(uint8_t); + *bsp_type = BSP_UINT8; + *element_size = sizeof(uint8_t); break; case mxINT64_CLASS: - bsp_type = BSP_INT64; - element_size = sizeof(int64_t); + *bsp_type = BSP_INT64; + *element_size = sizeof(int64_t); break; case mxINT32_CLASS: - bsp_type = BSP_INT32; - element_size = sizeof(int32_t); + *bsp_type = BSP_INT32; + *element_size = sizeof(int32_t); break; case mxINT16_CLASS: - bsp_type = BSP_INT16; - element_size = sizeof(int16_t); + *bsp_type = BSP_INT16; + *element_size = sizeof(int16_t); break; case mxINT8_CLASS: - bsp_type = BSP_INT8; - element_size = sizeof(int8_t); + *bsp_type = BSP_INT8; + *element_size = sizeof(int8_t); break; default: - return BSP_INVALID_TYPE; + return BSP_ERROR_TYPE; } } - if (is_empty) { + return BSP_SUCCESS; +} + +static inline bsp_error_t +matlab_to_bsp_array_allocator(const mxArray* mx_array, bsp_array_t* array, + bsp_allocator_t allocator) { + size_t size = mxGetNumberOfElements(mx_array); + bsp_type_t bsp_type; + size_t element_size; + bsp_error_t error = + matlab_bsp_array_type(mx_array, &bsp_type, &element_size); + if (error != BSP_SUCCESS) { + return error; + } + mxClassID class_id = mxGetClassID(mx_array); + bool is_complex = mxIsComplex(mx_array); + + if (mxIsEmpty(mx_array)) { array->data = NULL; array->size = 0; array->type = bsp_type; @@ -148,8 +177,7 @@ matlab_to_bsp_array_allocator(const mxArray* mx_array, bsp_array_t* array, return BSP_SUCCESS; } - bsp_error_t error = - bsp_construct_array_t_allocator(array, size, bsp_type, allocator); + error = bsp_construct_array_t_allocator(array, size, bsp_type, allocator); if (error != BSP_SUCCESS) { return error; } @@ -179,56 +207,45 @@ matlab_to_bsp_array_allocator(const mxArray* mx_array, bsp_array_t* array, 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; - } - +// Construct a non-owning Binsparse view of a real MATLAB numeric array. +// Complex MATLAB arrays use separate real and imaginary buffers with the MEX +// API used by these bindings, so they retain the existing owned conversion to +// C complex storage rather than exposing an invalid view. +static inline bsp_error_t matlab_to_bsp_array_view(const mxArray* mx_array, + bsp_array_t* array) { + bsp_type_t bsp_type; + size_t element_size; bsp_error_t error = - matlab_to_bsp_array_allocator(values_field, &matrix->values, allocator); + matlab_bsp_array_type(mx_array, &bsp_type, &element_size); + (void) element_size; 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; + if (mxIsComplex(mx_array)) { + return matlab_to_bsp_array_allocator(mx_array, array, + bsp_matlab_allocator); } - 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; + array->data = mxIsEmpty(mx_array) ? NULL : mxGetData(mx_array); + array->size = mxGetNumberOfElements(mx_array); + array->type = bsp_type; + array->allocator = bsp_matlab_view_allocator; + if (array->size > 0 && array->data == NULL) { + return BSP_ERROR_MEMORY; } + return BSP_SUCCESS; +} +static inline bsp_error_t +matlab_struct_to_bsp_matrix_metadata(const mxArray* mx_struct, + bsp_matrix_t* matrix) { 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; } @@ -239,13 +256,11 @@ static inline bsp_error_t matlab_struct_to_bsp_matrix_allocator( 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; } @@ -253,19 +268,16 @@ static inline bsp_error_t matlab_struct_to_bsp_matrix_allocator( 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; } @@ -279,6 +291,105 @@ static inline bsp_error_t matlab_struct_to_bsp_matrix_allocator( 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) { + if (!mx_struct || !mxIsStruct(mx_struct) || + mxGetNumberOfElements(mx_struct) != 1) { + return BSP_INVALID_STRUCTURE; + } + 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; + } + + error = matlab_struct_to_bsp_matrix_metadata(mx_struct, matrix); + if (error != BSP_SUCCESS) { + bsp_destroy_matrix_t(matrix); + } + return error; +} + +// Construct a matrix whose real arrays borrow the storage of the MATLAB BSP +// struct. The view is valid only while the MEX call's input remains alive and +// must be treated as read-only. Each array records its own ownership, allowing +// a minimizer to replace selected views with newly allocated MATLAB-owned +// arrays and then destroy the mixed-ownership matrix safely. +static inline bsp_error_t +matlab_struct_to_bsp_matrix_view(const mxArray* mx_struct, + bsp_matrix_t* matrix) { + if (!mx_struct || !mxIsStruct(mx_struct) || + mxGetNumberOfElements(mx_struct) != 1) { + return BSP_INVALID_STRUCTURE; + } + bsp_construct_default_matrix_t_allocator(matrix, + bsp_matlab_view_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) { + return BSP_INVALID_STRUCTURE; + } + + bsp_error_t error = matlab_to_bsp_array_view(values_field, &matrix->values); + if (error == BSP_SUCCESS) { + error = matlab_to_bsp_array_view(indices_0_field, &matrix->indices_0); + } + if (error == BSP_SUCCESS) { + error = matlab_to_bsp_array_view(indices_1_field, &matrix->indices_1); + } + if (error == BSP_SUCCESS) { + error = + matlab_to_bsp_array_view(pointers_to_1_field, &matrix->pointers_to_1); + } + if (error == BSP_SUCCESS) { + error = matlab_struct_to_bsp_matrix_metadata(mx_struct, matrix); + } + if (error != BSP_SUCCESS) { + bsp_destroy_matrix_t(matrix); + } + return error; +} + static inline bsp_error_t bsp_matrix_copy_with_allocator(const bsp_matrix_t* input, bsp_matrix_t* output, bsp_allocator_t allocator) { diff --git a/bindings/matlab/test_binsparse_from_ssmc.m b/bindings/matlab/test_binsparse_from_ssmc.m index 3d86e20..4939e1a 100644 --- a/bindings/matlab/test_binsparse_from_ssmc.m +++ b/bindings/matlab/test_binsparse_from_ssmc.m @@ -28,6 +28,18 @@ function test_binsparse_from_ssmc() error('Expected %d explicit zero values, got %d', nnz(Zeros), numel(zero_values)); end +assert(isa(mat.indices_1, 'uint8')); +assert(isa(mat.pointers_to_1, 'uint8')); + +% ISO values are detected before full conversion arrays are allocated. +iso = binsparse_from_ssmc(sparse([1 3], [2 4], [7 7], n, n), 'COO'); +assert(iso.is_iso && isequal(iso.values, 7)); +assert(isa(iso.indices_0, 'uint8') && isa(iso.indices_1, 'uint8')); + +% Explicit zeros preserve ISO only when every stored value is also zero. +zero_iso = binsparse_from_ssmc(sparse(n, n), Zeros, 'CSC'); +assert(zero_iso.is_iso && isequal(zero_iso.values, 0)); + fprintf('Test passed.\n'); end diff --git a/bindings/matlab/test_binsparse_matrix_views.m b/bindings/matlab/test_binsparse_matrix_views.m new file mode 100644 index 0000000..394f938 --- /dev/null +++ b/bindings/matlab/test_binsparse_matrix_views.m @@ -0,0 +1,98 @@ +function test_binsparse_matrix_views +%TEST_BINSPARSE_MATRIX_VIEWS exercise non-owning MATLAB BSP matrix views + +% SPDX-FileCopyrightText: 2026 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +required = {'binsparse_minimize_types', 'binsparse_write', 'binsparse_read'}; +for k = 1:numel(required) + if (exist(required{k}, 'file') ~= 3) + error('%s MEX function not found', required{k}); + end +end + +filename = [tempname '.bsp.h5']; +cleanup = onCleanup(@() delete_if_exists(filename)); %#ok + +% All four real arrays are borrowed. Values and indices are replaced during +% minimization, while the uint8 pointers remain a view until the output struct +% is constructed. Neither operation may alter or take ownership of the input. +input = make_matrix([1; 2; 4], uint64([]), uint64([0; 2; 1]), ... + uint8([0; 1; 2; 3]), 3, 3, 'CSC'); +output = binsparse_minimize_types(input); +assert(isa(input.values, 'double') && isequal(input.values, [1; 2; 4])); +assert(isa(input.indices_1, 'uint64') && ... + isequal(input.indices_1, uint64([0; 2; 1]))); +assert(isa(input.pointers_to_1, 'uint8') && ... + isequal(input.pointers_to_1, uint8([0; 1; 2; 3]))); +assert(isa(output.values, 'single')); +assert(isa(output.indices_1, 'uint8')); +assert(isa(output.pointers_to_1, 'uint8')); + +binsparse_write(filename, output); +assert(isequal(output.values, single([1; 2; 4]))); +assert(isequal(output.indices_1, uint8([0; 2; 1]))); +roundtrip = binsparse_read(filename); +assert(isequal(roundtrip.values, output.values)); +assert(isequal(roundtrip.indices_1, output.indices_1)); + +% If minimization changes nothing, every output data array is constructed from +% a view. The output must remain independently owned after the input is cleared. +unchanged = make_matrix(single([1; 2; 4]), uint8([0; 1; 2]), ... + uint8([0; 1; 2]), uint8([]), 3, 3, 'COO'); +unchanged_output = binsparse_minimize_types(unchanged); +clear unchanged +assert(isequal(unchanged_output.values, single([1; 2; 4]))); +assert(isequal(unchanged_output.indices_0, uint8([0; 1; 2]))); +assert(isequal(unchanged_output.indices_1, uint8([0; 1; 2]))); + +% Complex MATLAB storage is deliberately copied rather than viewed. Exercise +% both replacement by the minimizer and cleanup by the writer. +complex_input = make_matrix([1+2i; 3+4i], uint8([0; 1]), ... + uint8([0; 1]), uint8([]), 2, 2, 'COO'); +complex_output = binsparse_minimize_types(complex_input); +assert(isa(complex_input.values, 'double')); +assert(isequal(complex_input.values, [1+2i; 3+4i])); +assert(isa(complex_output.values, 'single')); +delete_if_exists(filename); +binsparse_write(filename, complex_output); +complex_roundtrip = binsparse_read(filename); +assert(isequal(complex_roundtrip.values, complex_output.values)); + +% Sparse arrays are invalid inside a raw BSP struct. Rejecting them prevents a +% view from treating a sparse value buffer as a dense numel-sized allocation. +invalid = make_matrix(sparse(eye(2)), uint8([0; 1]), ... + uint8([0; 1]), uint8([]), 2, 2, 'COO'); +assert_throws(@() binsparse_minimize_types(invalid), ... + 'BinSparse:ConversionError'); +assert_throws(@() binsparse_write(filename, invalid), ... + 'BinSparse:ConversionError'); + +fprintf('test_binsparse_matrix_views: all tests passed\n'); +end + +function matrix = make_matrix(values, indices_0, indices_1, pointers, ... + nrows, ncols, format) +matrix = struct('values', values, 'indices_0', indices_0, ... + 'indices_1', indices_1, 'pointers_to_1', pointers, ... + 'nrows', nrows, 'ncols', ncols, 'nnz', numel(values), ... + 'is_iso', false, 'format', format, 'structure', 'general'); +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', identifier); +end + +function delete_if_exists(filename) +if (exist(filename, 'file') == 2) + delete(filename); +end +end diff --git a/bindings/matlab/test_binsparse_write_ssmc_coo.m b/bindings/matlab/test_binsparse_write_ssmc_coo.m new file mode 100644 index 0000000..d8fe455 --- /dev/null +++ b/bindings/matlab/test_binsparse_write_ssmc_coo.m @@ -0,0 +1,93 @@ +function test_binsparse_write_ssmc_coo +%TEST_BINSPARSE_WRITE_SSMC_COO test the direct row-sorted COO writer + +% SPDX-FileCopyrightText: 2026 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +required = {'binsparse_write_ssmc_coo', 'binsparse_read'}; +for k = 1:numel(required) + if (exist(required{k}, 'file') ~= 3) + error('%s MEX function not found', required{k}); + end +end + +filename = [tempname '.bsp.h5']; +cleanup = onCleanup(@() delete_if_exists(filename)); %#ok + +A = sparse([3 1 2 1], [1 2 2 4], [5 6 7 8], 4, 4); +Zeros = sparse([4 2], [1 3], [1 1], 4, 4); +binsparse_write_ssmc_coo(filename, A, Zeros, '', ... + '{"metadata":{"role":"A"}}', 0); +raw = binsparse_read(filename); +assert(strcmp(raw.format, 'COO')); +assert(~raw.is_iso && isa(raw.values, 'single')); +assert(isa(raw.indices_0, 'uint8') && isa(raw.indices_1, 'uint8')); +assert(issortedrows([double(raw.indices_0), double(raw.indices_1)])); +assert(isequal(raw_to_sparse(raw), A)); +assert(nnz(raw.values == 0) == nnz(Zeros)); +assert(~isempty(regexp(h5readatt(filename, '/', 'binsparse'), ... + '"role"\s*:\s*"A"', 'once'))); + +% An ISO matrix stores one value even though all COO coordinates are retained. +delete_if_exists(filename); +iso = sparse([3 1 2], [1 2 4], [9 9 9], 4, 4); +binsparse_write_ssmc_coo(filename, iso, [], '', '{}', 9); +raw_iso = binsparse_read(filename); +assert(raw_iso.is_iso && numel(raw_iso.values) == 1); +assert(isequal(raw_to_sparse(raw_iso), iso)); + +% Dimensions, rather than the number of entries, select coordinate widths. +delete_if_exists(filename); +wide = sparse(70000, 70000, 1, 70000, 70000); +binsparse_write_ssmc_coo(filename, wide, [], '', '{}', 0); +raw_wide = binsparse_read(filename); +assert(isa(raw_wide.indices_0, 'uint32')); +assert(isa(raw_wide.indices_1, 'uint32')); + +% Complex values retain exactness and use the complex fallback safely. +delete_if_exists(filename); +complex_A = sparse([2 1], [1 2], [1+2i pi+4i], 2, 2); +binsparse_write_ssmc_coo(filename, complex_A, [], '', '{}', 0); +raw_complex = binsparse_read(filename); +assert(~raw_complex.is_iso && isa(raw_complex.values, 'double')); +assert(isequal(raw_to_sparse(raw_complex), complex_A)); + +% A named group appends to an existing primary file. +group_A = sparse([1 2], [2 1], [3 4], 2, 2); +binsparse_write_ssmc_coo(filename, group_A, [], 'aux', '{}', 0); +assert(isequal(raw_to_sparse(binsparse_read(filename, 'aux')), group_A)); + +duplicate = sparse(1, 1, 1, 2, 2); +assert_throws(@() binsparse_write_ssmc_coo( ... + filename, duplicate, duplicate, 'duplicate', '{}', 0), ... + 'BinSparse:DuplicateIndex'); + +fprintf('test_binsparse_write_ssmc_coo: all tests passed\n'); +end + +function A = raw_to_sparse(raw) +values = raw.values(:); +if (raw.is_iso) + values = repmat(values, raw.nnz, 1); +end +A = sparse(double(raw.indices_0) + 1, double(raw.indices_1) + 1, ... + double(values), raw.nrows, raw.ncols); +end + +function assert_throws(action, identifier) +try + action(); +catch me + assert(strcmp(me.identifier, identifier), ... + 'Expected %s, received %s', identifier, me.identifier); + return +end +error('Expected error %s', identifier); +end + +function delete_if_exists(filename) +if (exist(filename, 'file') == 2) + delete(filename); +end +end diff --git a/bindings/matlab/test_binsparse_write_ssmc_problem.m b/bindings/matlab/test_binsparse_write_ssmc_problem.m index 1f97af0..a96b275 100644 --- a/bindings/matlab/test_binsparse_write_ssmc_problem.m +++ b/bindings/matlab/test_binsparse_write_ssmc_problem.m @@ -9,6 +9,7 @@ function test_binsparse_write_ssmc_problem() required = {'binsparse_from_ssmc', 'binsparse_minimize_types', ... 'binsparse_write', 'binsparse_read', ... + 'binsparse_write_ssmc_coo', ... 'binsparse_write_string_dataset', 'binsparse_write_ssmc_problem'}; for i = 1:numel(required) if exist(required{i}, 'file') ~= 3 && exist(required{i}, 'file') ~= 2 @@ -68,6 +69,9 @@ function test_binsparse_write_ssmc_problem() % Read primary primary_bsp = binsparse_read(out_file); +assert(strcmp(primary_bsp.format, 'COO')); +assert(issortedrows([double(primary_bsp.indices_0), ... + double(primary_bsp.indices_1)])); primary_mat = bsp_to_matlab(primary_bsp); expected_primary = full(Problem.A); assert(matrices_equal(primary_mat, expected_primary), 'Primary matrix mismatch'); @@ -81,6 +85,7 @@ function test_binsparse_write_ssmc_problem() check_dense_group(out_file, 'D', Problem.aux.D); aux_sparse = binsparse_read(out_file, 'S'); +assert(strcmp(aux_sparse.format, 'COO')); 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'); diff --git a/include/binsparse/binsparse_cJSON.h b/include/binsparse/binsparse_cJSON.h index c74b9c4..e72b065 100644 --- a/include/binsparse/binsparse_cJSON.h +++ b/include/binsparse/binsparse_cJSON.h @@ -19,6 +19,8 @@ extern "C" { #endif +char* bsp_generate_json(bsp_matrix_t matrix, cJSON* user_json); + bsp_error_t bsp_write_matrix_cjson(const char* fname, bsp_matrix_t matrix, const char* group, cJSON* user_json, int compression_level);