diff --git a/.gitignore b/.gitignore index 77ff7ad..6b89fac 100644 --- a/.gitignore +++ b/.gitignore @@ -4,7 +4,6 @@ scripts venv -build/ build-*/ compile_flags.txt ._* diff --git a/CMakeLists.txt b/CMakeLists.txt index 5fd39e7..91247b2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,13 +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) @@ -26,6 +31,7 @@ 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) @@ -46,7 +52,8 @@ FetchContent_MakeAvailable(cJSON) 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} PRIVATE 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} @@ -57,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 @@ -64,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 @@ -76,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) @@ -99,6 +124,8 @@ if(CMAKE_PROJECT_NAME STREQUAL PROJECT_NAME) 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) diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..85e5a9b --- /dev/null +++ b/Makefile @@ -0,0 +1,71 @@ +#------------------------------------------------------------------------------- +# binsparse-reference-c/Makefile +#------------------------------------------------------------------------------- + +# SPDX-FileCopyrightText: 2024 Binsparse Developers +# +# SPDX-License-Identifier: BSD-3-Clause + +#------------------------------------------------------------------------------- + +# simple Makefile for binsparse-reference-c, relies on cmake to do the actual +# build. Use the CMAKE_OPTIONS argument to this Makefile to pass options to +# cmake. For example, to compile with 40 threads, use: +# +# make JOBS=40 + +JOBS ?= 8 + +default: library + +library: + ( cd build && cmake $(F) $(CMAKE_OPTIONS) .. && cmake --build . --config Release -j${JOBS} ) + +# install only in SuiteSparse/lib and SuiteSparse/include +local: + ( cd build && cmake $(F) $(CMAKE_OPTIONS) -USUITESPARSE_PKGFILEDIR -DSUITESPARSE_LOCAL_INSTALL=1 .. && cmake --build . --config Release -j${JOBS} ) + +# install only in /usr/local (default) +global: + ( cd build && cmake $(F) $(CMAKE_OPTIONS) -USUITESPARSE_PKGFILEDIR -DSUITESPARSE_LOCAL_INSTALL=0 .. && cmake --build . --config Release -j${JOBS} ) + +# compile with -g +debug: + ( cd build && cmake -DCMAKE_BUILD_TYPE=Debug $(F) $(CMAKE_OPTIONS) .. && cmake --build . --config Debug -j$(JOBS) ) + +# run the demos +demos: all + FIXME + +# just do 'make' in build; do not rerun the cmake script +remake: + ( cd build && cmake --build . -j$(JOBS) ) + +# just run cmake; do not compile +setup: + ( cd build && cmake $(F) $(CMAKE_OPTIONS) .. ) + +# build the static library +static: + ( cd build && cmake $(F) $(CMAKE_OPTIONS) -DBUILD_STATIC_LIBS=ON -DBUILD_SHARED_LIBS=OFF .. && cmake --build . --config Release -j$(JOBS) ) + +# installs to the install location defined by cmake, usually +# /usr/local/lib and /usr/local/include +install: + ( cd build && cmake --install . ) + +# create the doc +docs: + ( cd Doc && $(MAKE) ) + +# remove any installed libraries and #include files +uninstall: + - xargs rm < build/install_manifest.txt + +clean: distclean + +purge: distclean + +# remove all files not in the distribution +distclean: + - rm -rf build/* diff --git a/bindings/matlab/Contents.m b/bindings/matlab/Contents.m new file mode 100644 index 0000000..4cac8d2 --- /dev/null +++ b/bindings/matlab/Contents.m @@ -0,0 +1,31 @@ +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +% binsparse +% +% Files +% binsparse_read - read a sparse matrix from a binsparse hd5 file +% binsparse_write - write a matrix to a file in binsparse hd5 format +% binsparse_from_ssmc - convert SSMC A+Zeros to a Binsparse matrix struct +% binsparse_minimize_types - minimize value/index types in a Binsparse struct +% generate_bsp_from_ssmc - write SSMC problem to a Binsparse file +% +% bsp_matrix_create - SPDX-FileCopyrightText: 2024 Binsparse Developers +% bsp_matrix_info - SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% build_matlab_bindings - SPDX-FileCopyrightText: 2024 Binsparse Developers +% build_octave_bindings - SPDX-FileCopyrightText: 2024 Binsparse Developers +% compile_binsparse_read - SPDX-FileCopyrightText: 2024 Binsparse Developers +% compile_binsparse_read_octave - SPDX-FileCopyrightText: 2024 Binsparse Developers +% compile_binsparse_write - SPDX-FileCopyrightText: 2024 Binsparse Developers +% compile_binsparse_write_octave - SPDX-FileCopyrightText: 2024 Binsparse Developers +% compile_write_binsparse_from_matlab - SPDX-FileCopyrightText: 2024 Binsparse Developers +% compile_write_binsparse_from_matlab_octave - SPDX-FileCopyrightText: 2024 Binsparse Developers +% test_binsparse_read - SPDX-FileCopyrightText: 2024 Binsparse Developers +% test_binsparse_write - SPDX-FileCopyrightText: 2024 Binsparse Developers +% test_bsp_matrix_struct - SPDX-FileCopyrightText: 2024 Binsparse Developers +% test_binsparse_from_ssmc - SPDX-FileCopyrightText: 2024 Binsparse Developers +% test_binsparse_minimize_roundtrip - SPDX-FileCopyrightText: 2024 Binsparse Developers +% test_generate_bsp_from_ssmc - SPDX-FileCopyrightText: 2024 Binsparse Developers +% test_write_binsparse_from_matlab - SPDX-FileCopyrightText: 2024 Binsparse Developers diff --git a/bindings/matlab/README.md b/bindings/matlab/README.md index ff776f4..601fa11 100644 --- a/bindings/matlab/README.md +++ b/bindings/matlab/README.md @@ -57,6 +57,7 @@ mkoctfile --version ```matlab test_binsparse_read() test_binsparse_write() + test_write_binsparse_from_matlab() ``` #### Option 2: Octave (from within Octave) @@ -75,6 +76,7 @@ mkoctfile --version ```octave test_binsparse_read() test_binsparse_write() + test_write_binsparse_from_matlab() ``` #### Option 3: Octave (from command line) @@ -93,6 +95,7 @@ mkoctfile --version ```bash octave --eval "test_binsparse_read()" octave --eval "test_binsparse_write()" + octave --eval "test_write_binsparse_from_matlab()" ``` ## Usage Examples @@ -143,6 +146,25 @@ 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 optional parameters +write_binsparse_from_matlab(Problem, 'output.bsp.h5', 'my_group'); +write_binsparse_from_matlab(Problem, 'output.bsp.h5', 'my_group', '{"test": "metadata"}'); +write_binsparse_from_matlab(Problem, 'output.bsp.h5', 'my_group', '{"test": "metadata"}', 6); +``` + ### Error Handling The MEX functions include proper error handling: @@ -161,15 +183,19 @@ end |------|-------------| | `binsparse_read.c` | MEX function for reading Binsparse matrix files | | `binsparse_write.c` | MEX function for writing Binsparse matrix files | +| `write_binsparse_from_matlab.c` | MEX function for writing from SuiteSparse Matrix Collection format | | `build_matlab_bindings.m` | Main build script for MATLAB MEX functions | | `build_octave_bindings.m` | Main build script for Octave MEX functions | | `compile_binsparse_read.m` | Simple compilation script for read function (MATLAB) | | `compile_binsparse_write.m` | Simple compilation script for write function (MATLAB) | +| `compile_write_binsparse_from_matlab.m` | Simple compilation script for SuiteSparse write function (MATLAB) | | `compile_binsparse_read_octave.m` | Simple compilation script for read function (Octave) | | `compile_binsparse_write_octave.m` | Simple compilation script for write function (Octave) | +| `compile_write_binsparse_from_matlab_octave.m` | Simple compilation script for SuiteSparse write function (Octave) | | `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_write_binsparse_from_matlab.m` | Test script for SuiteSparse write functionality | | `bsp_matrix_create.m` | Utility function for creating matrix structs | | `bsp_matrix_info.m` | Utility function for displaying matrix information | | `README.md` | This documentation file | diff --git a/bindings/matlab/binsparse_from_ssmc.c b/bindings/matlab/binsparse_from_ssmc.c new file mode 100644 index 0000000..a8fd179 --- /dev/null +++ b/bindings/matlab/binsparse_from_ssmc.c @@ -0,0 +1,361 @@ +/* + * 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 "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 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, BSP_FLOAT64, + 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; + double* out_values = (double*) out->values.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])) { + out_values[out_ptr] = a->values[a_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])) { + out_values[out_ptr] = 0.0; + 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_DVEC) { + mexErrMsgIdAndTxt("BinSparse:InvalidFormat", + "Supported formats: CSC, CSR, COO, DMAT, 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, BSP_FLOAT64, + 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; + double* out_values = (double*) out->values.data; + + for (size_t i = 0; i < out->nnz; i++) { + out_values[i] = a->values[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) { + mexErrMsgIdAndTxt("BinSparse:InvalidFormat", + "Dense matrix requires DMAT 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 (mxIsComplex(mx_a) || !mxIsDouble(mx_a)) { + mexErrMsgIdAndTxt("BinSparse:InvalidMatrix", + "A must be a real 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))) { + mexErrMsgIdAndTxt("BinSparse:InvalidZeros", + "Zeros must be a real sparse double 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..45b0563 --- /dev/null +++ b/bindings/matlab/binsparse_from_ssmc.m @@ -0,0 +1,14 @@ +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +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, Zeros) +% matrix = binsparse_from_ssmc(A, Zeros, format) +% +% This is a thin wrapper over 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..e9f8d57 --- /dev/null +++ b/bindings/matlab/binsparse_minimize_types.c @@ -0,0 +1,523 @@ +/* + * 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 + +#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 bool float_values_are_int64(const double* values, size_t count) { + for (size_t i = 0; i < count; i++) { + double value = values[i]; + if (!isfinite(value)) { + return false; + } + if (value < (double) INT64_MIN || value > (double) INT64_MAX) { + return false; + } + int64_t ivalue = (int64_t) value; + if ((double) ivalue != value) { + return false; + } + } + return true; +} + +static bool float32_values_are_int64(const float* values, size_t count) { + for (size_t i = 0; i < count; i++) { + float value = values[i]; + if (!isfinite(value)) { + return false; + } + if (value < (float) INT64_MIN || value > (float) INT64_MAX) { + return false; + } + int64_t ivalue = (int64_t) value; + if ((float) ivalue != value) { + return false; + } + } + return true; +} + +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; + + if (float_values_are_int64(values, matrix->values.size)) { + bsp_array_t new_values; + bsp_error_t error = bsp_construct_array_t_allocator( + &new_values, matrix->values.size, BSP_INT64, bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + return; + } + + int64_t* n_values = (int64_t*) new_values.data; + for (size_t i = 0; i < matrix->values.size; i++) { + n_values[i] = (int64_t) values[i]; + } + + bsp_destroy_array_t(&matrix->values); + matrix->values = new_values; + minimize_int64_values(matrix); + return; + } + + 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_FLOAT32) { + float* values = (float*) matrix->values.data; + if (float32_values_are_int64(values, matrix->values.size)) { + bsp_array_t new_values; + bsp_error_t error = bsp_construct_array_t_allocator( + &new_values, matrix->values.size, BSP_INT64, bsp_matlab_allocator); + if (error != BSP_SUCCESS) { + return; + } + + int64_t* n_values = (int64_t*) new_values.data; + for (size_t i = 0; i < matrix->values.size; i++) { + n_values[i] = (int64_t) values[i]; + } + + bsp_destroy_array_t(&matrix->values); + matrix->values = new_values; + minimize_int64_values(matrix); + return; + } + } 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 index cf471fe..c9ae81d 100644 --- a/bindings/matlab/binsparse_read.c +++ b/bindings/matlab/binsparse_read.c @@ -19,14 +19,8 @@ #include #include -static inline void* bsp_matlab_malloc(size_t size) { - void* ptr = mxMalloc(size); - mexMakeMemoryPersistent(ptr); - return ptr; -} - -static const bsp_allocator_t bsp_matlab_allocator = { - .malloc = bsp_matlab_malloc, .free = mxFree}; +static const bsp_allocator_t bsp_matlab_allocator = {.malloc = mxMalloc, + .free = mxFree}; static inline mxClassID get_mxClassID(bsp_type_t type) { switch (type) { @@ -71,19 +65,26 @@ static inline mxComplexity get_mxComplexity(bsp_type_t type) { mxArray* bsp_array_to_matlab(bsp_array_t* array) { if (array->data == NULL || array->size == 0) { - // Return empty array - return mxCreateDoubleMatrix(1, 1, mxREAL); + // Return empty array of the correct class + mxClassID class_id = get_mxClassID(array->type); + if (class_id == mxUNKNOWN_CLASS) { + class_id = mxDOUBLE_CLASS; + } + return mxCreateNumericMatrix(0, 1, class_id, get_mxComplexity(array->type)); } if (get_mxClassID(array->type) == mxUNKNOWN_CLASS) { mexWarnMsgIdAndTxt("BinSparse:UnsupportedType", "Unsupported array type %d, returning empty array", (int) array->type); - return mxCreateDoubleMatrix(1, 1, mxREAL); + return mxCreateNumericMatrix(0, 1, mxDOUBLE_CLASS, mxREAL); } mxArray* mx_array = NULL; + // FIXME: do not use seperate real/imag arrays. Use the new MATLAB interleaved + // complex API. + if ((array->allocator.malloc == bsp_matlab_allocator.malloc && array->allocator.free == bsp_matlab_allocator.free) && get_mxComplexity(array->type) == mxREAL) { diff --git a/bindings/matlab/binsparse_read.m b/bindings/matlab/binsparse_read.m new file mode 100644 index 0000000..4dda8ce --- /dev/null +++ b/bindings/matlab/binsparse_read.m @@ -0,0 +1,18 @@ +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +function matrix = binsparse_read (filename, group) +%BINSPARSE_READ read a sparse matrix from a binsparse hd5 file +% +% Usage: +% matrix = binsparse_read (filename, group) +% +% FIXME: add documentation here +% +% Example: +% FIXME + +% FIXME add copyright + +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 index 9e56bd4..2057c14 100644 --- a/bindings/matlab/binsparse_write.c +++ b/bindings/matlab/binsparse_write.c @@ -21,217 +21,7 @@ #include #include -/** - * Convert MATLAB array to bsp_array_t - */ -bsp_error_t matlab_to_bsp_array(const mxArray* mx_array, bsp_array_t* array) { - if (mxIsEmpty(mx_array)) { - bsp_construct_default_array_t(array); - return BSP_SUCCESS; - } - - size_t size = mxGetNumberOfElements(mx_array); - mxClassID class_id = mxGetClassID(mx_array); - bool is_complex = mxIsComplex(mx_array); - - // Determine BSP type from MATLAB class - 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; - } - } - - // Allocate BSP array - bsp_error_t error = bsp_construct_array_t(array, size, bsp_type); - if (error != BSP_SUCCESS) { - return error; - } - - // Copy data - if (is_complex) { - // Handle complex numbers: interleave real/imaginary parts - if (class_id == mxDOUBLE_CLASS) { - double* real_data = mxGetPr(mx_array); - double* imag_data = mxGetPi(mx_array); - double* out_data = (double*) array->data; - for (size_t i = 0; i < size; i++) { - out_data[2 * i] = real_data[i]; // Real part - out_data[2 * i + 1] = imag_data[i]; // Imaginary part - } - } else { // mxSINGLE_CLASS - float* real_data = (float*) mxGetData(mx_array); - float* imag_data = (float*) mxGetImagData(mx_array); - float* out_data = (float*) array->data; - for (size_t i = 0; i < size; i++) { - out_data[2 * i] = real_data[i]; // Real part - out_data[2 * i + 1] = imag_data[i]; // Imaginary part - } - } - } else { - // Simple copy for real types - memcpy(array->data, mxGetData(mx_array), size * element_size); - } - - return BSP_SUCCESS; -} - -/** - * Convert MATLAB struct to bsp_matrix_t - */ -bsp_error_t matlab_struct_to_bsp_matrix(const mxArray* mx_struct, - bsp_matrix_t* matrix) { - bsp_construct_default_matrix_t(matrix); - - // Extract and convert arrays - 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; - error = matlab_to_bsp_array(values_field, &matrix->values); - if (error != BSP_SUCCESS) { - bsp_destroy_matrix_t(matrix); - return error; - } - - error = matlab_to_bsp_array(indices_0_field, &matrix->indices_0); - if (error != BSP_SUCCESS) { - bsp_destroy_matrix_t(matrix); - return error; - } - - error = matlab_to_bsp_array(indices_1_field, &matrix->indices_1); - if (error != BSP_SUCCESS) { - bsp_destroy_matrix_t(matrix); - return error; - } - - error = matlab_to_bsp_array(pointers_to_1_field, &matrix->pointers_to_1); - if (error != BSP_SUCCESS) { - bsp_destroy_matrix_t(matrix); - return error; - } - - // Extract scalar fields - 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); - - // Extract format string - 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; - } - - // Extract structure string - 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; // Default fallback - } - - return BSP_SUCCESS; -} +#include "matlab_bsp_helpers.h" /** * Main MEX function entry point @@ -274,7 +64,8 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { } // Convert MATLAB struct to bsp_matrix_t - error = matlab_struct_to_bsp_matrix(prhs[1], &matrix); + error = matlab_struct_to_bsp_matrix_allocator(prhs[1], &matrix, + bsp_matlab_allocator); if (error != BSP_SUCCESS) { mxFree(filename); mexErrMsgIdAndTxt("BinSparse:ConversionError", diff --git a/bindings/matlab/binsparse_write.m b/bindings/matlab/binsparse_write.m new file mode 100644 index 0000000..9d4542d --- /dev/null +++ b/bindings/matlab/binsparse_write.m @@ -0,0 +1,16 @@ +function binsparse_write(filename, matrix, group, json_string, compression_level) +%BINSPARSE_WRITE write a matrix to a file in binsparse hd5 format +% +% Usage: +% binsparse_write (filename, matrix, group, json_string, compression_level) +% +% FIXME: add documentation here +% +% Example: +% +% FIXME + +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% SPDX-License-Identifier: BSD-3-Clause + +error ('binsparse_write mexFunction not found; compile with build_matlab_bindings first') ; diff --git a/bindings/matlab/build_matlab_bindings.m b/bindings/matlab/build_matlab_bindings.m index 8816d0f..61b4814 100644 --- a/bindings/matlab/build_matlab_bindings.m +++ b/bindings/matlab/build_matlab_bindings.m @@ -55,7 +55,8 @@ function build_matlab_bindings(varargin) 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\n'); +fprintf(' test_binsparse_write()\n'); +fprintf(' test_write_binsparse_from_matlab()\n\n'); end @@ -68,7 +69,8 @@ function build_matlab_bindings(varargin) if success fprintf('MEX compiler found: %s\n', cc(1).Name); end - catch + catch me + me success = false; end end @@ -94,7 +96,9 @@ function compile_mex_functions(paths, verbose) % Compile all MEX functions % List of MEX functions to compile - mex_files = {'binsparse_read.c', 'binsparse_write.c'}; + mex_files = {'binsparse_read.c', 'binsparse_write.c', ... + 'binsparse_from_ssmc.c', 'binsparse_minimize_types.c', ... + 'write_binsparse_from_matlab.c'}; fprintf('Compiling MEX functions...\n'); @@ -108,21 +112,36 @@ function compile_mex_functions(paths, verbose) fprintf(' Compiling %s... ', mex_file); % Prepare MEX command with library linking + % FIXME: use .so not .a lib_dir = fullfile(paths.binsparse_root, 'build'); - lib_path = fullfile(lib_dir, 'libbinsparse.a'); - cjson_lib = fullfile(lib_dir, '_deps', 'cjson-build', 'libcjson.so'); - - mex_args = {'-I', paths.include_dir, mex_file, lib_path, cjson_lib, '-lhdf5_serial'}; - if verbose - mex_args = [mex_args, {'-v'}]; + 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=' ; + end + rpath = sprintf (' -Wl,%s''''%s'''' ', rpath, lib_dir) ; + rpath = [' LDFLAGS=''$LDFLAGS -fPIC ' rpath ' '' '] ; + +% mex_args = {'-I', paths.include_dir, mex_file, lib_path, cjson_lib, '-lhdf5_serial'}; + paths.include_dir + mex_args = sprintf ('mex -g -I%s %s %s %s %s -lhdf5_serial', ... + paths.include_dir, rpath, lib_path, mex_file, cjson_lib) ; +% if verbose + if 0 +% mex_args = [mex_args, {'-v'}]; + mex_args = [mex_args, ' -v'] ; end try - mex(mex_args{:}); + mex_args + eval (mex_args) ; fprintf('SUCCESS\n'); - catch ME + catch me + me fprintf('FAILED\n'); - fprintf(' Error: %s\n', ME.message); + fprintf(' Error: %s\n', me.message); end end end diff --git a/bindings/matlab/build_octave_bindings.m b/bindings/matlab/build_octave_bindings.m index 87ba30f..4ab99dd 100644 --- a/bindings/matlab/build_octave_bindings.m +++ b/bindings/matlab/build_octave_bindings.m @@ -60,7 +60,8 @@ function build_octave_bindings(varargin) 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\n'); +fprintf(' test_binsparse_write()\n'); +fprintf(' test_write_binsparse_from_matlab()\n\n'); end @@ -103,7 +104,9 @@ 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'}; + mex_files = {'binsparse_read.c', 'binsparse_write.c', ... + 'binsparse_from_ssmc.c', 'binsparse_minimize_types.c', ... + 'write_binsparse_from_matlab.c'}; fprintf('Compiling MEX functions with mkoctfile...\n'); diff --git a/bindings/matlab/compile_binsparse_write.m b/bindings/matlab/compile_binsparse_write.m index d62bb08..6eb86b9 100644 --- a/bindings/matlab/compile_binsparse_write.m +++ b/bindings/matlab/compile_binsparse_write.m @@ -34,7 +34,8 @@ function compile_binsparse_write() try % Compile with linking - mex('-I', include_dir, 'binsparse_write.c', lib_path, cjson_lib, '-lhdf5_serial', '-v'); + mex('-I', include_dir, 'binsparse_write.c', lib_path, cjson_lib, ... + '-lhdf5_serial', '-v'); fprintf('Successfully compiled binsparse_write!\n'); % Test if it loads diff --git a/bindings/matlab/compile_octave.sh b/bindings/matlab/compile_octave.sh index 641edf2..527def2 100755 --- a/bindings/matlab/compile_octave.sh +++ b/bindings/matlab/compile_octave.sh @@ -137,7 +137,7 @@ if [ "$CLEAN" = true ]; then fi # List of MEX files to compile -MEX_FILES=("binsparse_read.c" "binsparse_write.c") +MEX_FILES=("binsparse_read.c" "binsparse_write.c" "binsparse_from_ssmc.c" "binsparse_minimize_types.c" "write_binsparse_from_matlab.c") print_info "Compiling MEX functions..." diff --git a/bindings/matlab/compile_write_binsparse_from_matlab.m b/bindings/matlab/compile_write_binsparse_from_matlab.m new file mode 100644 index 0000000..2e17bcb --- /dev/null +++ b/bindings/matlab/compile_write_binsparse_from_matlab.m @@ -0,0 +1,173 @@ +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +function compile_write_binsparse_from_matlab(varargin) +% COMPILE_WRITE_BINSPARSE_FROM_MATLAB - Compile the write_binsparse_from_matlab MEX function +% +% This script compiles the write_binsparse_from_matlab MEX function for MATLAB. +% It automatically detects include paths and links against the Binsparse library. +% +% Usage: +% compile_write_binsparse_from_matlab() % Standard compilation +% compile_write_binsparse_from_matlab('verbose') % Verbose compilation output +% compile_write_binsparse_from_matlab('debug') % Debug build with symbols +% +% Prerequisites: +% - MATLAB with working MEX compiler (run 'mex -setup' if needed) +% - Binsparse C library headers (in ../../include/) +% - Compiled Binsparse library (in ../../build/) + +% Parse input arguments +verbose = any(strcmpi(varargin, 'verbose')); +debug_build = any(strcmpi(varargin, 'debug')); + +fprintf('=== Compiling write_binsparse_from_matlab MEX function ===\n\n'); + +% Check if source file exists +source_file = 'write_binsparse_from_matlab.c'; +if ~exist(source_file, 'file') + error('Source file not found: %s\nEnsure you are in the correct directory.', source_file); +end + +% Check MEX compiler +if ~check_mex_compiler() + error('MEX compiler not properly configured. Run "mex -setup" first.'); +end + +% Get build 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); + fprintf(' Library dir: %s\n', paths.lib_dir); + fprintf(' Root dir: %s\n', paths.binsparse_root); + fprintf('\n'); +end + +% Compile the MEX function +compile_mex_function(source_file, paths, verbose, debug_build); + +fprintf('\n=== Compilation Complete ===\n'); +fprintf('Test the function with:\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.lib_dir = fullfile(paths.binsparse_root, 'build'); + + 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 + + if ~exist(paths.lib_dir, 'dir') + error('Binsparse build directory not found: %s\nEnsure the library has been built first.', paths.lib_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 + + % Check for compiled library + lib_path = fullfile(paths.lib_dir, 'libbinsparse.a'); + if ~exist(lib_path, 'file') + error('Binsparse library not found: %s\nEnsure the library has been built first.', lib_path); + end +end + +function compile_mex_function(source_file, paths, verbose, debug_build) + % Compile the MEX function with appropriate flags and libraries + + fprintf('Compiling %s... ', source_file); + + % Prepare MEX command with library linking + lib_path = fullfile(paths.lib_dir, 'libbinsparse.a'); + cjson_lib = fullfile(paths.lib_dir, '_deps', 'cjson-build', 'libcjson.so'); + + % Check if we're on macOS and adjust cjson library path + if ismac + % On macOS, cjson might have a different extension or location + cjson_alternatives = { + fullfile(paths.lib_dir, '_deps', 'cjson-build', 'libcjson.dylib'), + fullfile(paths.lib_dir, '_deps', 'cjson-build', 'libcjson.a') + }; + + for i = 1:length(cjson_alternatives) + if exist(cjson_alternatives{i}, 'file') + cjson_lib = cjson_alternatives{i}; + break; + end + end + end + + % Build MEX arguments + mex_args = {'-I', paths.include_dir, source_file, lib_path}; + + % Add cjson library if it exists + if exist(cjson_lib, 'file') + mex_args{end+1} = cjson_lib; + else + warning('cjson library not found at expected location: %s', cjson_lib); + fprintf(' Continuing without cjson library...\n'); + end + + % Add HDF5 library + mex_args{end+1} = '-lhdf5_serial'; + + % Add optional flags + if verbose + mex_args{end+1} = '-v'; + end + + if debug_build + mex_args = [mex_args, {'-g', '-DDEBUG'}]; + fprintf('(debug build) '); + end + + if verbose + fprintf('\n MEX command: '); + for i = 1:length(mex_args) + fprintf('%s ', mex_args{i}); + end + fprintf('\n'); + end + + try + mex(mex_args{:}); + fprintf('SUCCESS\n'); + catch ME + fprintf('FAILED\n'); + fprintf(' Error: %s\n', ME.message); + + % Provide troubleshooting suggestions + fprintf('\nTroubleshooting suggestions:\n'); + fprintf(' 1. Ensure the Binsparse library has been built: cd ../../build && make\n'); + fprintf(' 2. Check that MEX compiler is configured: mex -setup\n'); + fprintf(' 3. Verify HDF5 development libraries are installed\n'); + fprintf(' 4. Try running with verbose flag for more details\n'); + + rethrow(ME); + end +end diff --git a/bindings/matlab/compile_write_binsparse_from_matlab_octave.m b/bindings/matlab/compile_write_binsparse_from_matlab_octave.m new file mode 100644 index 0000000..8bf6003 --- /dev/null +++ b/bindings/matlab/compile_write_binsparse_from_matlab_octave.m @@ -0,0 +1,150 @@ +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +function compile_write_binsparse_from_matlab_octave(varargin) +% COMPILE_WRITE_BINSPARSE_FROM_MATLAB_OCTAVE - Compile the write_binsparse_from_matlab MEX function for Octave +% +% This script compiles the write_binsparse_from_matlab MEX function for Octave using mkoctfile. +% It automatically detects include paths and links against the Binsparse library. +% +% Usage: +% compile_write_binsparse_from_matlab_octave() % Standard compilation +% compile_write_binsparse_from_matlab_octave('verbose') % Verbose compilation output +% +% Prerequisites: +% - GNU Octave with mkoctfile +% - Binsparse C library headers (in ../../include/) +% - Compiled Binsparse library (in ../../build/) + +% Parse input arguments +verbose = any(strcmpi(varargin, 'verbose')); + +fprintf('=== Compiling write_binsparse_from_matlab MEX function for Octave ===\n\n'); + +% Check if we're running in Octave +if ~is_octave() + warning('This script is designed for Octave. For MATLAB, use compile_write_binsparse_from_matlab.m'); +end + +% Check if source file exists +source_file = 'write_binsparse_from_matlab.c'; +if ~exist(source_file, 'file') + error('Source file not found: %s\nEnsure you are in the correct directory.', source_file); +end + +% Check mkoctfile availability +if ~check_mkoctfile() + error('mkoctfile not found. Please ensure Octave is properly installed.'); +end + +% Get build 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); + fprintf(' Library dir: %s\n', paths.lib_dir); + fprintf(' Root dir: %s\n', paths.binsparse_root); + fprintf('\n'); +end + +% Compile the MEX function +compile_octave_function(source_file, paths, verbose); + +fprintf('\n=== Compilation Complete ===\n'); +fprintf('Test the function with:\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.lib_dir = fullfile(paths.binsparse_root, 'build'); + + 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 + + if ~exist(paths.lib_dir, 'dir') + error('Binsparse build directory not found: %s\nEnsure the library has been built first.', paths.lib_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 + + % Check for compiled library + lib_path = fullfile(paths.lib_dir, 'libbinsparse.a'); + if ~exist(lib_path, 'file') + error('Binsparse library not found: %s\nEnsure the library has been built first.', lib_path); + end +end + +function compile_octave_function(source_file, paths, verbose) + % Compile the MEX function using mkoctfile + + fprintf('Compiling %s with mkoctfile... ', source_file); + + % Prepare mkoctfile command with library linking + include_flag = sprintf('-I%s', paths.include_dir); + lib_path = fullfile(paths.lib_dir, 'libbinsparse.a'); + cjson_lib_dir = fullfile(paths.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, source_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, source_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); + + % Provide troubleshooting suggestions + fprintf('\nTroubleshooting suggestions:\n'); + fprintf(' 1. Ensure the Binsparse library has been built: cd ../../build && make\n'); + fprintf(' 2. Check that Octave development packages are installed\n'); + fprintf(' 3. Verify HDF5 development libraries are installed\n'); + fprintf(' 4. Try running with verbose flag for more details\n'); + + error('Compilation failed'); + end +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..b5e17a6 --- /dev/null +++ b/bindings/matlab/generate_bsp_from_ssmc.m @@ -0,0 +1,108 @@ +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +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 + +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); +binsparse_write(output_filename, mat, '', [], 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 + +% Metadata handling not yet implemented +fprintf('Note: metadata handling is not implemented yet.\n'); + +end + +function handle_aux_entry(name, value, output_filename, format, compression_level) + if ischar(value) || (isstring(value) && isscalar(value)) + fprintf('Note: aux entry "%s" is a string and is currently ignored.\n', name); + return; + 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, [], compression_level); +end + +function fmt = dense_format_for(value) + if isvector(value) + fmt = 'DVEC'; + else + fmt = 'DMAT'; + end +end diff --git a/bindings/matlab/matlab_bsp_helpers.h b/bindings/matlab/matlab_bsp_helpers.h new file mode 100644 index 0000000..18a7871 --- /dev/null +++ b/bindings/matlab/matlab_bsp_helpers.h @@ -0,0 +1,475 @@ +/* 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 + +typedef struct { + double* values; + mwIndex* rowind; + mwIndex* colptr; + size_t nrows; + size_t ncols; + size_t nnz; +} 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 (mxIsComplex(mx_matrix)) { + return -1; + } + + if (!mxIsDouble(mx_matrix)) { + return -1; + } + + csc_matrix->nrows = mxGetM(mx_matrix); + csc_matrix->ncols = mxGetN(mx_matrix); + csc_matrix->nnz = mxGetNzmax(mx_matrix); + + csc_matrix->values = mxGetPr(mx_matrix); + csc_matrix->rowind = mxGetIr(mx_matrix); + csc_matrix->colptr = mxGetJc(mx_matrix); + + if (!csc_matrix->values || !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; + } + + 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* out_data = (double*) array->data; + for (size_t i = 0; i < size; i++) { + out_data[2 * i] = real_data[i]; + out_data[2 * i + 1] = imag_data[i]; + } + } else { + float* real_data = (float*) mxGetData(mx_array); + float* imag_data = (float*) mxGetImagData(mx_array); + float* out_data = (float*) array->data; + for (size_t i = 0; i < size; i++) { + out_data[2 * i] = real_data[i]; + out_data[2 * i + 1] = imag_data[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, get_mxClassID(array->type), get_mxComplexity(array->type)); + 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* in_data = (float*) 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] = in_data[2 * i]; + imag_data[i] = in_data[2 * i + 1]; + } + } else { + double* in_data = (double*) 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] = in_data[2 * i]; + imag_data[i] = in_data[2 * i + 1]; + } + } + } + } + + 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..9dfe57b --- /dev/null +++ b/bindings/matlab/test_binsparse_from_ssmc.m @@ -0,0 +1,33 @@ +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +function test_binsparse_from_ssmc() +% TEST_BINSPARSE_FROM_SSMC - Basic test for binsparse_from_ssmc MEX function + +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..1ba5a07 --- /dev/null +++ b/bindings/matlab/test_binsparse_minimize_roundtrip.m @@ -0,0 +1,67 @@ +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +function test_binsparse_minimize_roundtrip() +% TEST_BINSPARSE_MINIMIZE_ROUNDTRIP - Test SSMC conversion + type minimization + +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_roundtrip_dir.m b/bindings/matlab/test_binsparse_roundtrip_dir.m new file mode 100644 index 0000000..a7dad81 --- /dev/null +++ b/bindings/matlab/test_binsparse_roundtrip_dir.m @@ -0,0 +1,156 @@ +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +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. + +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_generate_bsp_from_ssmc.m b/bindings/matlab/test_generate_bsp_from_ssmc.m new file mode 100644 index 0000000..78a50cf --- /dev/null +++ b/bindings/matlab/test_generate_bsp_from_ssmc.m @@ -0,0 +1,145 @@ +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +function test_generate_bsp_from_ssmc() +% TEST_GENERATE_BSP_FROM_SSMC - End-to-end test for generate_bsp_from_ssmc + +fprintf('=== Testing generate_bsp_from_ssmc ===\n\n'); + +required = {'binsparse_from_ssmc', 'binsparse_minimize_types', ... + 'binsparse_write', 'binsparse_read', '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.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 = 'ignored.txt'; + +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'); + +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 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.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..ad8a865 --- /dev/null +++ b/bindings/matlab/test_write_binsparse_from_matlab.m @@ -0,0 +1,149 @@ +% SPDX-FileCopyrightText: 2024 Binsparse Developers +% +% SPDX-License-Identifier: BSD-3-Clause + +function test_write_binsparse_from_matlab() +% TEST_WRITE_BINSPARSE_FROM_MATLAB - Test the write_binsparse_from_matlab MEX function +% +% This function tests the write_binsparse_from_matlab MEX function by creating +% sample SuiteSparse Matrix Collection problem structs and attempting to write +% them to Binsparse format files. + +fprintf('=== Testing write_binsparse_from_matlab MEX function ===\n\n'); + +try + % Test 1: Check if MEX function exists and is callable + fprintf('Test 1: Checking MEX function availability\n'); + if exist('write_binsparse_from_matlab', 'file') ~= 3 + error('write_binsparse_from_matlab MEX function not found. Please compile it first.'); + end + fprintf('✓ MEX function found\n\n'); + + % Test 2: Test with minimal arguments (should fail) + fprintf('Test 2: Testing error handling with insufficient arguments\n'); + try + write_binsparse_from_matlab(); + error('Expected error for insufficient arguments'); + catch ME + if ~isempty(strfind(ME.identifier, 'BinSparse:InvalidArgs')) + fprintf('✓ Correctly handled insufficient arguments\n'); + else + fprintf('✗ Unexpected error: %s\n', ME.message); + end + end + fprintf('\n'); + + % Test 3: Test with invalid problem struct + fprintf('Test 3: Testing error handling with invalid problem struct\n'); + try + write_binsparse_from_matlab(42, 'test.bsp.h5'); + error('Expected error for invalid problem struct'); + catch ME + if ~isempty(strfind(ME.identifier, 'BinSparse:InvalidProblemStruct')) + fprintf('✓ Correctly handled invalid problem struct\n'); + else + fprintf('✗ Unexpected error: %s\n', ME.message); + end + end + fprintf('\n'); + + % Test 4: Create a basic SuiteSparse Matrix Collection problem struct + fprintf('Test 4: Testing with basic SuiteSparse problem struct\n'); + problem = create_basic_problem_struct(); + problem + + fprintf('Testing skeleton implementation with basic problem struct:\n'); + write_binsparse_from_matlab(problem, 'test_basic.bsp.h5'); + fprintf('✓ Basic test completed successfully\n\n'); + + % Test 5: Test with all optional arguments + fprintf('Test 5: Testing with all optional arguments\n'); + problem_extended = create_extended_problem_struct(); + + fprintf('Testing with all arguments:\n'); + write_binsparse_from_matlab(problem_extended, 'test_extended.bsp.h5', ... + 'my_group', '{"test": "metadata"}', 6); + fprintf('✓ Extended test completed successfully\n\n'); + + % Test 6: Test with real sparse matrix + fprintf('Test 6: Testing with real sparse matrix\n'); + problem_sparse = create_sparse_problem_struct(); + + fprintf('Testing with sparse matrix:\n'); + write_binsparse_from_matlab(problem_sparse, 'test_sparse.bsp.h5'); + fprintf('✓ Sparse matrix test completed successfully\n\n'); + + fprintf('=== All tests passed! ===\n'); + fprintf('Note: This is testing the skeleton implementation only.\n'); + fprintf('The actual file writing functionality needs to be implemented.\n\n'); + +catch ME + fprintf('✗ Test failed: %s\n', ME.message); + fprintf('Stack trace:\n'); + for i = 1:length(ME.stack) + fprintf(' %s (line %d)\n', ME.stack(i).name, ME.stack(i).line); + end +end + +end + +function problem = create_basic_problem_struct() + % Create a minimal SuiteSparse Matrix Collection problem struct + problem = struct(); + problem.name = 'test_basic_matrix'; + problem.A = speye(3); % 3x3 identity matrix + problem.title = 'Basic Test Matrix'; + problem.kind = 'test matrix'; +end + +function problem = create_extended_problem_struct() + % Create an extended SuiteSparse Matrix Collection problem struct + problem = struct(); + problem.name = 'test_extended_matrix'; + problem.title = 'Extended Test Matrix with Metadata'; + problem.kind = 'artificial/test'; + problem.A = sparse([1 2 3], [1 2 3], [1.5 2.5 3.5], 4, 4); + problem.notes = 'This is a test matrix for validation'; + problem.author = 'Test Suite'; + problem.date = datestr(now); + problem.editor = 'MATLAB'; + + % Add some additional fields that might be present + problem.id = 12345; + problem.group = 'Test'; + problem.num_rows = size(problem.A, 1); + problem.num_cols = size(problem.A, 2); + problem.nnz = nnz(problem.A); + problem.pattern_symmetry = 1.0; + problem.numerical_symmetry = 1.0; + problem.type = 'real'; + problem.structure = 'unsymmetric'; + problem.sprank = rank(full(problem.A)); +end + +function problem = create_sparse_problem_struct() + % Create a problem struct with a more complex sparse matrix + problem = struct(); + problem.name = 'test_sparse_matrix'; + problem.title = 'Sparse Test Matrix'; + problem.kind = 'test/sparse'; + + % Create a 10x10 sparse matrix with some pattern + n = 10; + [i, j] = meshgrid(1:n, 1:n); + mask = abs(i - j) <= 2; % Pentadiagonal pattern + rows = i(mask); + cols = j(mask); + vals = randn(length(rows), 1); % Random values + + problem.A = sparse(rows, cols, vals, n, n); + problem.notes = sprintf('Random %dx%d pentadiagonal matrix with %d non-zeros', ... + n, n, nnz(problem.A)); + + % Add matrix properties: FIXME: why? + problem.num_rows = n; + problem.num_cols = n; + problem.nnz = nnz(problem.A); + problem.type = 'real'; + problem.structure = 'unsymmetric'; +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..8623575 --- /dev/null +++ b/bindings/matlab/write_binsparse_from_matlab.c @@ -0,0 +1,467 @@ +/* + * SPDX-FileCopyrightText: 2024 Binsparse Developers + * + * SPDX-License-Identifier: BSD-3-Clause + */ + +/** + * write_binsparse_from_matlab.c - Write SuiteSparse Matrix Collection Problem + * struct to Binsparse format + * + * This MEX function takes a SuiteSparse Matrix Collection problem struct and + * converts it to the Binsparse format, writing the result to a specified file. + * + * Usage in MATLAB/Octave: + * write_binsparse_from_matlab(problem_struct, filename) + * write_binsparse_from_matlab(problem_struct, filename, group) + * write_binsparse_from_matlab(problem_struct, filename, group, json_metadata) + * write_binsparse_from_matlab(problem_struct, filename, group, json_metadata, + * compression_level) + * + * Arguments: + * problem_struct - SuiteSparse Matrix Collection problem struct with fields: + * .name - Problem name (string) + * .A - Sparse matrix (MATLAB sparse matrix) + * .Zeros - Sparse matrix with pattern of explicit zeros in the problem + * .title - Problem title (optional) + * .kind - Problem kind (optional) + * .notes - Additional notes (optional) + * filename - Output filename for the Binsparse file + * group - Optional HDF5 group name (default: 'default') + * json_metadata - Optional JSON metadata string + * compression_level - Optional compression level (0-9, default: 1) + */ + +#include "mex.h" +#include +#include +#include +#include + +static const bsp_allocator_t bsp_matlab_allocator = {.malloc = mxMalloc, + .free = mxFree}; + +typedef struct { + double* values; + mwIndex* rowind; + mwIndex* colptr; + size_t nrows; + size_t ncols; + size_t nnz; +} matlab_csc_t; + +int extract_matlab_csc(const mxArray* mx_matrix, matlab_csc_t* csc_matrix) { + // Validate input + if (!mx_matrix) { + mexPrintf("Error: NULL matrix pointer\n"); + return -1; + } + + if (!mxIsSparse(mx_matrix)) { + mexPrintf("Error: Matrix is not sparse\n"); + return -1; + } + + if (mxIsComplex(mx_matrix)) { + mexPrintf("Error: Complex matrices not yet supported\n"); + return -1; + } + + // Extract matrix dimensions + csc_matrix->nrows = mxGetM(mx_matrix); + csc_matrix->ncols = mxGetN(mx_matrix); + csc_matrix->nnz = mxGetNzmax(mx_matrix); + + // Get pointers to MATLAB's internal CSC data + // Note: MATLAB stores sparse matrices in CSC format internally + csc_matrix->values = mxGetPr(mx_matrix); // Non-zero values + csc_matrix->rowind = mxGetIr(mx_matrix); // Row indices (0-based) + csc_matrix->colptr = mxGetJc(mx_matrix); // Column pointers + + // Validate that we got valid pointers + if (!csc_matrix->values || !csc_matrix->rowind || !csc_matrix->colptr) { + mexPrintf("Error: Failed to extract CSC data from MATLAB matrix\n"); + return -1; + } + + // The actual number of non-zeros might be less than nzmax + if (csc_matrix->ncols > 0) { + csc_matrix->nnz = + csc_matrix->colptr[csc_matrix->ncols]; // Last element of colptr gives + // actual nnz + } + + mexPrintf("Extracted CSC matrix: %zu x %zu with %zu non-zeros\n", + csc_matrix->nrows, csc_matrix->ncols, csc_matrix->nnz); + + return 0; // Success +} + +bsp_matrix_t merge_csc_with_zeros(matlab_csc_t matrix, matlab_csc_t zeros) { + + mexPrintf("Merging %zu x %zu matrix (%zu nnz) with zeros pattern (%zu nnz)\n", + matrix.nrows, matrix.ncols, matrix.nnz, zeros.nnz); + + // If there are no zeros, just construct matrix based on `matrix` and return. + if (zeros.nnz == 0) { + bsp_matrix_t result; + bsp_construct_default_matrix_t_allocator(&result, bsp_matlab_allocator); + + result.nrows = matrix.nrows; + result.ncols = matrix.ncols; + result.nnz = matrix.nnz; + result.format = BSP_CSC; + result.structure = BSP_GENERAL; + result.is_iso = false; + + bsp_construct_array_t_allocator(&result.values, matrix.nnz, BSP_FLOAT64, + bsp_matlab_allocator); + bsp_construct_array_t_allocator(&result.pointers_to_1, matrix.ncols + 1, + BSP_UINT64, bsp_matlab_allocator); + bsp_construct_array_t_allocator(&result.indices_1, matrix.nnz, BSP_UINT64, + bsp_matlab_allocator); + + double* result_values = (double*) result.values.data; + /* + OLD: + mwIndex* result_rowind = (mwIndex*) result.pointers_to_1.data; + mwIndex* result_colptr = (mwIndex*) result.indices_1.data; + */ + // NEW: + mwIndex* result_rowind = (mwIndex*) result.indices_1.data; + mwIndex* result_colptr = (mwIndex*) result.pointers_to_1.data; + + for (size_t i = 0; i < result.values.size; i++) { + result_values[i] = matrix.values[i]; + } + + for (size_t i = 0; i < result.pointers_to_1.size; i++) { + result_colptr[i] = matrix.colptr[i]; + } + + for (size_t i = 0; i < result.indices_1.size; i++) { + result_rowind[i] = matrix.rowind[i]; + } + + return result; + } + + // Initialize result, which will hold the merged `matrix` and `zeros`. + bsp_matrix_t result; + bsp_construct_default_matrix_t_allocator(&result, bsp_matlab_allocator); + + // Set up result matrix metadata + result.nrows = matrix.nrows; + result.ncols = matrix.ncols; + result.nnz = matrix.nnz + zeros.nnz; + result.format = BSP_CSC; + result.structure = BSP_GENERAL; + result.is_iso = false; + + assert(sizeof(mwIndex) == sizeof(uint64_t)); + + bsp_construct_array_t_allocator(&result.values, matrix.nnz + zeros.nnz, + BSP_FLOAT64, bsp_matlab_allocator); + bsp_construct_array_t_allocator(&result.pointers_to_1, matrix.ncols + 1, + BSP_UINT64, bsp_matlab_allocator); + bsp_construct_array_t_allocator(&result.indices_1, matrix.nnz + zeros.nnz, + BSP_UINT64, bsp_matlab_allocator); + + double* result_values = (double*) result.values.data; + mwIndex* result_colptr = (mwIndex*) result.pointers_to_1.data; + mwIndex* result_rowind = (mwIndex*) result.indices_1.data; + + // Set colptr for result + result_colptr[0] = 0; + + for (mwIndex j = 1; j < matrix.ncols + 1; j++) { + mwIndex row_nnz_matrix = matrix.colptr[j] - matrix.colptr[j - 1]; + mwIndex row_nnz_zeros = zeros.colptr[j] - zeros.colptr[j - 1]; + result_colptr[j] = result_colptr[j - 1] + row_nnz_matrix + row_nnz_zeros; + } + + // FIXME: this produces a bsp_matrix_t with row indices out of order. + // Is that OK? + + for (mwIndex j = 0; j < result.ncols; j++) { + mwIndex result_i_ptr = result_colptr[j]; + + for (mwIndex matrix_i_ptr = matrix.colptr[j]; + matrix_i_ptr < matrix.colptr[j + 1]; matrix_i_ptr++) { + result_values[result_i_ptr] = matrix.values[matrix_i_ptr]; + result_rowind[result_i_ptr] = matrix.rowind[matrix_i_ptr]; + + result_i_ptr++; + } + + for (mwIndex zeros_i_ptr = zeros.colptr[j]; + zeros_i_ptr < zeros.colptr[j + 1]; zeros_i_ptr++) { + result_values[result_i_ptr] = + zeros.values[zeros_i_ptr]; // FIXME: should be 0 + result_rowind[result_i_ptr] = zeros.rowind[zeros_i_ptr]; + + result_i_ptr++; + } + + assert(result_i_ptr == result_colptr[j + 1]); + } + + mexPrintf("Successfully created merged bsp_matrix_t in CSC format\n"); + return result; +} + +/** + * Print information about a SuiteSparse Matrix Collection problem struct + */ +void print_problem_info(const mxArray* problem_struct) { + mexPrintf("=== SuiteSparse Matrix Collection Problem Information ===\n"); + + // Check if input is a struct + if (!mxIsStruct(problem_struct)) { + mexPrintf("Error: Input is not a struct\n"); + return; + } + + mexPrintf("Number of fields: %d\n", mxGetNumberOfFields(problem_struct)); + mexPrintf("Number of elements: %d\n", + (int) mxGetNumberOfElements(problem_struct)); + + // Print field names + mexPrintf("Field names:\n"); + for (int i = 0; i < mxGetNumberOfFields(problem_struct); i++) { + const char* field_name = mxGetFieldNameByNumber(problem_struct, i); + mexPrintf(" [%d] %s\n", i, field_name); + + // Get field value and print basic info + mxArray* field_value = mxGetFieldByNumber(problem_struct, 0, i); + if (field_value) { + if (mxIsChar(field_value)) { + char* str_value = mxArrayToString(field_value); + if (str_value) { + mexPrintf(" (string): \"%s\"\n", str_value); + mxFree(str_value); + } + } else if (mxIsSparse(field_value)) { + mexPrintf(" (sparse matrix): %dx%d with %d non-zeros\n", + (int) mxGetM(field_value), (int) mxGetN(field_value), + (int) mxGetNzmax(field_value)); + } else if (mxIsNumeric(field_value)) { + mexPrintf(" (numeric): %dx%d %s array\n", + (int) mxGetM(field_value), (int) mxGetN(field_value), + mxGetClassName(field_value)); + } else { + mexPrintf(" (%s): %dx%d\n", mxGetClassName(field_value), + (int) mxGetM(field_value), (int) mxGetN(field_value)); + } + } else { + mexPrintf(" (null)\n"); + } + } + mexPrintf("=========================================================\n\n"); +} + +/** + * Main MEX function entry point + */ +void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { + char* filename = NULL; + char* group = NULL; + char* json_metadata = NULL; + int compression_level = 0; // Default compression + + // Check input arguments + if (nrhs < 2 || nrhs > 5) { + mexErrMsgIdAndTxt( + "BinSparse:InvalidArgs", + "Usage: write_binsparse_from_matlab(problem_struct, filename [, group " + "[, json_metadata [, compression_level]]])"); + } + + if (nlhs > 0) { + mexErrMsgIdAndTxt("BinSparse:TooManyOutputs", + "No output arguments expected"); + } + + mexPrintf("Number of input arguments: %d\n", nrhs); + mexPrintf("Number of output arguments: %d\n", nlhs); + + // Validate and process problem struct (first argument) + if (!mxIsStruct(prhs[0])) { + mexErrMsgIdAndTxt("BinSparse:InvalidProblemStruct", + "First argument must be a SuiteSparse Matrix Collection " + "problem struct"); + } + + mxArray* mx_problem = mxGetField(prhs[0], 0, "Problem"); + + if ((mx_problem == NULL) || !mxIsStruct(mx_problem)) { + mexErrMsgIdAndTxt("BinSparse:InvalidProblemStruct", + "First argument must be a SuiteSparse Matrix Collection " + "problem struct"); + } + + // Extract sparse matrix from problem.A field + mxArray* mx_matrix = mxGetField(mx_problem, 0, "A"); + + if (!mx_matrix) { + mexErrMsgIdAndTxt("BinSparse:InvalidProblemStruct", + "First argument must be a SuiteSparse Matrix Collection " + "problem struct"); + } + + if (!mxIsSparse(mx_matrix)) { + mexErrMsgIdAndTxt( + "BinSparse:InvalidProblemStruct", + "Unable to extract Matlab sparse matrix --- not a matrix"); + } + + mexPrintf("Found sparse matrix in problem.A field\n"); + matlab_csc_t csc_matrix = {0}; + int rv = extract_matlab_csc(mx_matrix, &csc_matrix); + + if (rv != 0) { + mexErrMsgIdAndTxt("BinSparse:InvalidProblemStruct", + "Unable to extract Matlab sparse matrix"); + } + + mxArray* mx_zeros_matrix = mxGetField(mx_problem, 0, "Zeros"); + + matlab_csc_t zeros_csc_matrix = {0}; + + if (mx_zeros_matrix) { + if (!mxIsSparse(mx_zeros_matrix)) { + mexErrMsgIdAndTxt("BinSparse:InvalidProblemStruct", + "Zeros matrix exists but is not sparse"); + } + + int rv = extract_matlab_csc(mx_zeros_matrix, &zeros_csc_matrix); + + if (rv != 0) { + mexErrMsgIdAndTxt("BinSparse:InvalidProblemStruct", + "Unable to extract Zeros sparse matrix"); + } + } + + bsp_matrix_t merged_matrix = + merge_csc_with_zeros(csc_matrix, zeros_csc_matrix); + + mexPrintf("Merged matrix!\n"); + + mexPrintf("Merged matrix is %zu x %zu and has %zu nnz.\n", + merged_matrix.nrows, merged_matrix.ncols, merged_matrix.nnz); + + bsp_destroy_matrix_t(&merged_matrix); + + mexPrintf("\nAnalyzing problem struct:\n"); + print_problem_info(prhs[0]); + + // Get filename (second argument) + if (!mxIsChar(prhs[1])) { + mexErrMsgIdAndTxt("BinSparse:InvalidFilename", "Filename must be a string"); + } + + filename = mxArrayToString(prhs[1]); + if (!filename) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to convert filename string"); + } + mexPrintf("Output filename: \"%s\"\n", filename); + + // Get optional group name (third argument) + if (nrhs >= 3 && !mxIsEmpty(prhs[2])) { + if (!mxIsChar(prhs[2])) { + mexErrMsgIdAndTxt("BinSparse:InvalidGroup", + "Group name must be a string"); + } + + group = mxArrayToString(prhs[2]); + if (!group) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to convert group string"); + } + mexPrintf("HDF5 group: \"%s\"\n", group); + } else { + mexPrintf("HDF5 group: (default)\n"); + } + + // Get optional JSON metadata (fourth argument) + if (nrhs >= 4 && !mxIsEmpty(prhs[3])) { + if (!mxIsChar(prhs[3])) { + mexErrMsgIdAndTxt("BinSparse:InvalidJSON", + "JSON metadata must be a string"); + } + + json_metadata = mxArrayToString(prhs[3]); + if (!json_metadata) { + mexErrMsgIdAndTxt("BinSparse:MemoryError", + "Failed to convert JSON metadata string"); + } + mexPrintf("JSON metadata: \"%s\"\n", json_metadata); + } else { + mexPrintf("JSON metadata: (none)\n"); + } + + // Get optional compression level (fifth argument) + if (nrhs >= 5 && !mxIsEmpty(prhs[4])) { + if (!mxIsNumeric(prhs[4]) || mxIsComplex(prhs[4]) || + mxGetNumberOfElements(prhs[4]) != 1) { + mexErrMsgIdAndTxt("BinSparse:InvalidCompression", + "Compression level must be a scalar integer"); + } + + compression_level = (int) mxGetScalar(prhs[4]); + mexPrintf("Compression level: %d\n", compression_level); + } else { + mexPrintf("Compression level: %d (default)\n", compression_level); + } + + mexPrintf("\n=== IMPLEMENTATION STATUS ===\n"); + + // Extract sparse matrix from problem.A field + mxArray* matrix_field = mxGetField(prhs[0], 0, "A"); + if (matrix_field && mxIsSparse(matrix_field)) { + mexPrintf("✓ Found sparse matrix in problem.A field\n"); + + // Test the matlab_csc_to_bsp function + mexPrintf("Testing matlab_csc_to_bsp function:\n"); + matlab_csc_t csc; + extract_matlab_csc(matrix_field, &csc); + + if (csc.values && csc.rowind && csc.colptr) { + mexPrintf("✓ Successfully extracted CSC data:\n"); + mexPrintf(" - Dimensions: %zu x %zu\n", csc.nrows, csc.ncols); + mexPrintf(" - Non-zeros: %zu\n", csc.nnz); + mexPrintf(" - Values pointer: %p\n", (void*) csc.values); + mexPrintf(" - Row indices pointer: %p\n", (void*) csc.rowind); + mexPrintf(" - Column pointers: %p\n", (void*) csc.colptr); + + // Show first few values as example + if (csc.nnz > 0) { + mexPrintf(" - First few values: "); + size_t show_count = (csc.nnz < 5) ? csc.nnz : 5; + for (size_t i = 0; i < show_count; i++) { + mexPrintf("%.6g ", csc.values[i]); + } + if (csc.nnz > 5) + mexPrintf("..."); + mexPrintf("\n"); + } + } else { + mexPrintf("✗ Failed to extract CSC data\n"); + } + } else { + mexPrintf("✗ No sparse matrix found in problem.A field\n"); + } + + mexPrintf("\nTODO: Convert MATLAB sparse matrix to bsp_matrix_t\n"); + mexPrintf("TODO: Add metadata from problem struct to JSON\n"); + mexPrintf("TODO: Call bsp_write_matrix() to write the file\n"); + mexPrintf("=====================================\n\n"); + + mexPrintf("Function completed successfully (skeleton mode)\n"); + + // Clean up allocated strings + mxFree(json_metadata); + mxFree(group); + mxFree(filename); +} diff --git a/examples/benchmark_write.c b/examples/benchmark_write.c index d2e5a60..eaed5ab 100644 --- a/examples/benchmark_write.c +++ b/examples/benchmark_write.c @@ -5,6 +5,7 @@ */ #include +#include #include #include #include @@ -70,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) { diff --git a/examples/bsp2mtx.c b/examples/bsp2mtx.c index 8b40c32..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; } diff --git a/include/binsparse/array.h b/include/binsparse/array.h index b3b7db2..39604aa 100644 --- a/include/binsparse/array.h +++ b/include/binsparse/array.h @@ -25,6 +25,7 @@ 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; } diff --git a/include/binsparse/convert_matrix.h b/include/binsparse/convert_matrix.h index 93c96ef..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); @@ -21,7 +38,7 @@ static inline bsp_matrix_t bsp_convert_matrix(bsp_matrix_t matrix, if (matrix.format == BSP_CSR) { // Convert CSR -> COOR bsp_matrix_t result; - bsp_construct_default_matrix_t(&result); + bsp_construct_default_matrix_t_allocator(&result, allocator); result.format = BSP_COOR; @@ -38,11 +55,11 @@ static inline bsp_matrix_t bsp_convert_matrix(bsp_matrix_t matrix, bsp_type_t index_type = bsp_pick_integer_type(max_dim); - bsp_error_t error = - bsp_copy_construct_array_t(&result.values, 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(&empty_result); + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); return empty_result; } @@ -51,20 +68,21 @@ static inline bsp_matrix_t bsp_convert_matrix(bsp_matrix_t matrix, // we might upcast. if (index_type == matrix.indices_1.type) { - error = bsp_copy_construct_array_t(&result.indices_1, 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(&empty_result); + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); return empty_result; } } else { - error = bsp_construct_array_t(&result.indices_1, 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(&empty_result); + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); return empty_result; } @@ -75,12 +93,13 @@ static inline bsp_matrix_t bsp_convert_matrix(bsp_matrix_t matrix, } } - error = bsp_construct_array_t(&result.indices_0, 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(&empty_result); + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); return empty_result; } @@ -93,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); } @@ -102,8 +283,10 @@ 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_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 { @@ -111,7 +294,7 @@ static inline bsp_matrix_t bsp_convert_matrix(bsp_matrix_t matrix, // Convert COOR -> CSR bsp_matrix_t result; - bsp_construct_default_matrix_t(&result); + bsp_construct_default_matrix_t_allocator(&result, allocator); result.format = BSP_CSR; @@ -137,30 +320,30 @@ 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. - bsp_error_t error = - bsp_copy_construct_array_t(&result.values, 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(&empty_result); + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); return empty_result; } if (index_type == matrix.indices_1.type) { - error = - bsp_copy_construct_array_t(&result.indices_1, 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(&empty_result); + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); return empty_result; } } else { - error = - bsp_construct_array_t(&result.indices_1, 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(&empty_result); + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); return empty_result; } @@ -171,13 +354,13 @@ static inline bsp_matrix_t bsp_convert_matrix(bsp_matrix_t matrix, } } - error = bsp_construct_array_t(&result.pointers_to_1, 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(&empty_result); + bsp_construct_default_matrix_t_allocator(&empty_result, allocator); return empty_result; } @@ -203,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); @@ -210,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/hdf5_wrapper.h b/include/binsparse/hdf5_wrapper.h index 16ff7a7..f8b44df 100644 --- a/include/binsparse/hdf5_wrapper.h +++ b/include/binsparse/hdf5_wrapper.h @@ -293,6 +293,7 @@ bsp_read_attribute_allocator(char** string, hid_t f, const char* label, *string = (char*) allocator.malloc(size + 1); H5Aread(attribute, strtype, *string); + (*string)[size] = '\0'; H5Aclose(attribute); H5Tclose(strtype); diff --git a/include/binsparse/matrix.h b/include/binsparse/matrix.h index bedb3ed..bc54f85 100644 --- a/include/binsparse/matrix.h +++ b/include/binsparse/matrix.h @@ -30,7 +30,7 @@ typedef struct bsp_matrix_t { static inline void bsp_construct_default_matrix_t_allocator(bsp_matrix_t* matrix, bsp_allocator_t allocator) { - bsp_construct_default_array_t(&matrix->values); + 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); 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/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 e06a82b..d9fb788 100644 --- a/src/read_matrix.c +++ b/src/read_matrix.c @@ -154,12 +154,15 @@ bsp_error_t bsp_read_matrix_from_group_parallel(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) { + printf("construct default allocator:\n"); bsp_construct_default_matrix_t_allocator(matrix, allocator); char* json_string; + printf("read attr allocator:\n"); bsp_error_t error = bsp_read_attribute_allocator( &json_string, f, (char*) "binsparse", allocator); if (error != BSP_SUCCESS) { + printf("fail here!:\n"); return error; } @@ -222,7 +225,7 @@ bsp_error_t bsp_read_matrix_from_group_allocator(bsp_matrix_t* matrix, hid_t f, error = bsp_read_array_allocator(&matrix->values, f, (char*) "values", allocator); if (error != BSP_SUCCESS) { - free(json_string); + allocator.free(json_string); return error; } @@ -247,7 +250,7 @@ bsp_error_t bsp_read_matrix_from_group_allocator(bsp_matrix_t* matrix, hid_t f, error = bsp_read_array_allocator(&matrix->indices_0, f, (char*) "indices_0", allocator); if (error != BSP_SUCCESS) { - free(json_string); + allocator.free(json_string); bsp_destroy_array_t(&matrix->values); return error; } @@ -257,7 +260,7 @@ bsp_error_t bsp_read_matrix_from_group_allocator(bsp_matrix_t* matrix, hid_t f, error = bsp_read_array_allocator(&matrix->indices_1, f, (char*) "indices_1", allocator); if (error != BSP_SUCCESS) { - free(json_string); + allocator.free(json_string); bsp_destroy_array_t(&matrix->values); bsp_destroy_array_t(&matrix->indices_0); return error; @@ -268,7 +271,7 @@ bsp_error_t bsp_read_matrix_from_group_allocator(bsp_matrix_t* matrix, hid_t f, error = bsp_read_array_allocator(&matrix->pointers_to_1, f, (char*) "pointers_to_1", allocator); if (error != BSP_SUCCESS) { - free(json_string); + allocator.free(json_string); bsp_destroy_array_t(&matrix->values); bsp_destroy_array_t(&matrix->indices_0); bsp_destroy_array_t(&matrix->indices_1); @@ -284,7 +287,7 @@ bsp_error_t bsp_read_matrix_from_group_allocator(bsp_matrix_t* matrix, hid_t f, } cJSON_Delete(j); - free(json_string); + allocator.free(json_string); return BSP_SUCCESS; } diff --git a/src/write_matrix.c b/src/write_matrix.c index 7d1956d..ede7788 100644 --- a/src/write_matrix.c +++ b/src/write_matrix.c @@ -21,10 +21,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* item_copy = cJSON_Duplicate(item, 1); // 1 = deep copy - cJSON_AddItemToObject(j, item->string, item_copy); + 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); @@ -136,7 +138,13 @@ bsp_error_t bsp_write_matrix_to_group_cjson(hid_t f, bsp_matrix_t matrix, 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 = cJSON_Parse(user_json); + 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); @@ -181,7 +189,13 @@ bsp_error_t bsp_write_matrix_cjson(const char* fname, bsp_matrix_t matrix, 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 = cJSON_Parse(user_json); + 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);