Skip to content

Commit 4b4d3f0

Browse files
BenBrockBenBrock
committed
Begin implementing generate_bsp_from_ssmc, add support for
minimization.
1 parent 37203f6 commit 4b4d3f0

14 files changed

Lines changed: 1731 additions & 216 deletions

‎bindings/matlab/Contents.m‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,9 @@
33
% Files
44
% binsparse_read - read a sparse matrix from a binsparse hd5 file
55
% binsparse_write - write a matrix to a file in binsparse hd5 format
6+
% binsparse_from_ssmc - convert SSMC A+Zeros to a Binsparse matrix struct
7+
% binsparse_minimize_types - minimize value/index types in a Binsparse struct
8+
% generate_bsp_from_ssmc - write SSMC problem to a Binsparse file
69
%
710
% bsp_matrix_create - SPDX-FileCopyrightText: 2024 Binsparse Developers
811
% bsp_matrix_info - SPDX-FileCopyrightText: 2024 Binsparse Developers
@@ -18,4 +21,7 @@
1821
% test_binsparse_read - SPDX-FileCopyrightText: 2024 Binsparse Developers
1922
% test_binsparse_write - SPDX-FileCopyrightText: 2024 Binsparse Developers
2023
% test_bsp_matrix_struct - SPDX-FileCopyrightText: 2024 Binsparse Developers
24+
% test_binsparse_from_ssmc - SPDX-FileCopyrightText: 2024 Binsparse Developers
25+
% test_binsparse_minimize_roundtrip - SPDX-FileCopyrightText: 2024 Binsparse Developers
26+
% test_generate_bsp_from_ssmc - SPDX-FileCopyrightText: 2024 Binsparse Developers
2127
% test_write_binsparse_from_matlab - SPDX-FileCopyrightText: 2024 Binsparse Developers

‎bindings/matlab/binsparse_from_ssmc.c‎

Lines changed: 360 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,360 @@
1+
/*
2+
* SPDX-FileCopyrightText: 2024 Binsparse Developers
3+
*
4+
* SPDX-License-Identifier: BSD-3-Clause
5+
*/
6+
7+
/**
8+
* binsparse_from_ssmc.c - Convert SuiteSparse matrix A + Zeros to a Binsparse
9+
* MATLAB struct.
10+
*
11+
* Usage in MATLAB/Octave:
12+
* matrix = binsparse_from_ssmc(A)
13+
* matrix = binsparse_from_ssmc(A, Zeros)
14+
* matrix = binsparse_from_ssmc(A, Zeros, format)
15+
* matrix = binsparse_from_ssmc(A, format)
16+
*
17+
* Inputs:
18+
* A - sparse or dense matrix (MATLAB)
19+
* Zeros - optional sparse matrix representing explicit zeros (same size as
20+
* A) format - optional string: default 'COO' for sparse, 'DMAT'/'DVEC' for
21+
* dense
22+
*
23+
* Output:
24+
* MATLAB struct compatible with binsparse_write.
25+
*/
26+
27+
#include "mex.h"
28+
#include <binsparse/binsparse.h>
29+
#include <string.h>
30+
31+
#include "matlab_bsp_helpers.h"
32+
33+
static bool array_uses_allocator(bsp_array_t array, bsp_allocator_t allocator) {
34+
if (array.size == 0 || array.data == NULL) {
35+
return true;
36+
}
37+
return array.allocator.malloc == allocator.malloc &&
38+
array.allocator.free == allocator.free;
39+
}
40+
41+
static bool matrix_uses_allocator(const bsp_matrix_t* matrix,
42+
bsp_allocator_t allocator) {
43+
return array_uses_allocator(matrix->values, allocator) &&
44+
array_uses_allocator(matrix->indices_0, allocator) &&
45+
array_uses_allocator(matrix->indices_1, allocator) &&
46+
array_uses_allocator(matrix->pointers_to_1, allocator);
47+
}
48+
49+
static bsp_error_t construct_array_with_allocator(bsp_array_t* array,
50+
size_t size, bsp_type_t type,
51+
bsp_allocator_t allocator) {
52+
if (size == 0) {
53+
array->data = NULL;
54+
array->size = 0;
55+
array->type = type;
56+
array->allocator = allocator;
57+
return BSP_SUCCESS;
58+
}
59+
return bsp_construct_array_t_allocator(array, size, type, allocator);
60+
}
61+
62+
static void build_csc_merged(const matlab_csc_t* a, const matlab_csc_t* z,
63+
bsp_matrix_t* out) {
64+
bsp_error_t error;
65+
66+
bsp_construct_default_matrix_t_allocator(out, bsp_matlab_allocator);
67+
out->nrows = a->nrows;
68+
out->ncols = a->ncols;
69+
out->nnz = a->nnz + z->nnz;
70+
out->format = BSP_CSC;
71+
out->structure = BSP_GENERAL;
72+
out->is_iso = false;
73+
74+
error = construct_array_with_allocator(&out->values, out->nnz, BSP_FLOAT64,
75+
bsp_matlab_allocator);
76+
if (error != BSP_SUCCESS) {
77+
mexErrMsgIdAndTxt("BinSparse:MemoryError",
78+
"Failed to allocate values array");
79+
}
80+
81+
error = construct_array_with_allocator(&out->indices_1, out->nnz, BSP_UINT64,
82+
bsp_matlab_allocator);
83+
if (error != BSP_SUCCESS) {
84+
mexErrMsgIdAndTxt("BinSparse:MemoryError",
85+
"Failed to allocate indices array");
86+
}
87+
88+
error = construct_array_with_allocator(&out->pointers_to_1, out->ncols + 1,
89+
BSP_UINT64, bsp_matlab_allocator);
90+
if (error != BSP_SUCCESS) {
91+
mexErrMsgIdAndTxt("BinSparse:MemoryError",
92+
"Failed to allocate pointers array");
93+
}
94+
95+
uint64_t* out_colptr = (uint64_t*) out->pointers_to_1.data;
96+
uint64_t* out_rowind = (uint64_t*) out->indices_1.data;
97+
double* out_values = (double*) out->values.data;
98+
99+
out_colptr[0] = 0;
100+
for (mwIndex j = 0; j < a->ncols; j++) {
101+
mwIndex a_count = a->colptr[j + 1] - a->colptr[j];
102+
mwIndex z_count = z->colptr[j + 1] - z->colptr[j];
103+
out_colptr[j + 1] = out_colptr[j] + a_count + z_count;
104+
}
105+
106+
for (mwIndex j = 0; j < a->ncols; j++) {
107+
mwIndex a_ptr = a->colptr[j];
108+
mwIndex a_end = a->colptr[j + 1];
109+
mwIndex z_ptr = z->colptr[j];
110+
mwIndex z_end = z->colptr[j + 1];
111+
uint64_t out_ptr = out_colptr[j];
112+
113+
while (a_ptr < a_end || z_ptr < z_end) {
114+
if (z_ptr >= z_end ||
115+
(a_ptr < a_end && a->rowind[a_ptr] < z->rowind[z_ptr])) {
116+
out_values[out_ptr] = a->values[a_ptr];
117+
out_rowind[out_ptr] = (uint64_t) a->rowind[a_ptr];
118+
a_ptr++;
119+
} else if (a_ptr >= a_end ||
120+
(z_ptr < z_end && z->rowind[z_ptr] < a->rowind[a_ptr])) {
121+
out_values[out_ptr] = 0.0;
122+
out_rowind[out_ptr] = (uint64_t) z->rowind[z_ptr];
123+
z_ptr++;
124+
} else {
125+
mexErrMsgIdAndTxt("BinSparse:DuplicateIndex",
126+
"Duplicate indices between A and Zeros");
127+
}
128+
out_ptr++;
129+
}
130+
131+
if (out_ptr != out_colptr[j + 1]) {
132+
mexErrMsgIdAndTxt("BinSparse:InternalError",
133+
"Merged column counts do not match");
134+
}
135+
}
136+
}
137+
138+
static bsp_matrix_format_t parse_format(int nrhs, const mxArray* prhs[]) {
139+
const mxArray* format_arg = NULL;
140+
if (nrhs >= 3 && !mxIsEmpty(prhs[2])) {
141+
format_arg = prhs[2];
142+
} else if (nrhs == 2 && mxIsChar(prhs[1])) {
143+
format_arg = prhs[1];
144+
}
145+
146+
if (!format_arg) {
147+
return BSP_COO;
148+
}
149+
150+
if (!mxIsChar(format_arg)) {
151+
mexErrMsgIdAndTxt("BinSparse:InvalidFormat", "Format must be a string");
152+
}
153+
154+
char* format_str = mxArrayToString(format_arg);
155+
if (!format_str) {
156+
mexErrMsgIdAndTxt("BinSparse:MemoryError", "Failed to read format string");
157+
}
158+
159+
bsp_matrix_format_t format = bsp_get_matrix_format(format_str);
160+
mxFree(format_str);
161+
162+
if (format != BSP_CSC && format != BSP_CSR && format != BSP_COO &&
163+
format != BSP_COOR && format != BSP_DMAT && format != BSP_DVEC) {
164+
mexErrMsgIdAndTxt("BinSparse:InvalidFormat",
165+
"Supported formats: CSC, CSR, COO, DMAT, DVEC");
166+
}
167+
168+
return format;
169+
}
170+
171+
static void build_csc_from_a(const matlab_csc_t* a, bsp_matrix_t* out) {
172+
bsp_error_t error;
173+
174+
bsp_construct_default_matrix_t_allocator(out, bsp_matlab_allocator);
175+
out->nrows = a->nrows;
176+
out->ncols = a->ncols;
177+
out->nnz = a->nnz;
178+
out->format = BSP_CSC;
179+
out->structure = BSP_GENERAL;
180+
out->is_iso = false;
181+
182+
error = construct_array_with_allocator(&out->values, out->nnz, BSP_FLOAT64,
183+
bsp_matlab_allocator);
184+
if (error != BSP_SUCCESS) {
185+
mexErrMsgIdAndTxt("BinSparse:MemoryError",
186+
"Failed to allocate values array");
187+
}
188+
189+
error = construct_array_with_allocator(&out->indices_1, out->nnz, BSP_UINT64,
190+
bsp_matlab_allocator);
191+
if (error != BSP_SUCCESS) {
192+
mexErrMsgIdAndTxt("BinSparse:MemoryError",
193+
"Failed to allocate indices array");
194+
}
195+
196+
error = construct_array_with_allocator(&out->pointers_to_1, out->ncols + 1,
197+
BSP_UINT64, bsp_matlab_allocator);
198+
if (error != BSP_SUCCESS) {
199+
mexErrMsgIdAndTxt("BinSparse:MemoryError",
200+
"Failed to allocate pointers array");
201+
}
202+
203+
uint64_t* out_colptr = (uint64_t*) out->pointers_to_1.data;
204+
uint64_t* out_rowind = (uint64_t*) out->indices_1.data;
205+
double* out_values = (double*) out->values.data;
206+
207+
for (size_t i = 0; i < out->nnz; i++) {
208+
out_values[i] = a->values[i];
209+
out_rowind[i] = (uint64_t) a->rowind[i];
210+
}
211+
212+
for (size_t i = 0; i < out->ncols + 1; i++) {
213+
out_colptr[i] = (uint64_t) a->colptr[i];
214+
}
215+
}
216+
217+
static void build_dense_matrix(const mxArray* mx_a, bsp_matrix_t* out,
218+
bsp_matrix_format_t format) {
219+
if (!mxIsNumeric(mx_a)) {
220+
mexErrMsgIdAndTxt("BinSparse:InvalidMatrix", "Dense input must be numeric");
221+
}
222+
223+
bool is_vector = (mxGetM(mx_a) == 1) || (mxGetN(mx_a) == 1);
224+
225+
if (is_vector && format != BSP_DVEC) {
226+
mexErrMsgIdAndTxt("BinSparse:InvalidFormat",
227+
"Dense vector requires DVEC format");
228+
}
229+
230+
if (!is_vector && format != BSP_DMAT) {
231+
mexErrMsgIdAndTxt("BinSparse:InvalidFormat",
232+
"Dense matrix requires DMAT format");
233+
}
234+
235+
bsp_construct_default_matrix_t_allocator(out, bsp_matlab_allocator);
236+
out->format = format;
237+
out->structure = BSP_GENERAL;
238+
out->is_iso = false;
239+
240+
size_t nrows = mxGetM(mx_a);
241+
size_t ncols = mxGetN(mx_a);
242+
size_t total = mxGetNumberOfElements(mx_a);
243+
244+
if (is_vector) {
245+
out->nrows = total;
246+
out->ncols = 1;
247+
} else {
248+
out->nrows = nrows;
249+
out->ncols = ncols;
250+
}
251+
252+
out->nnz = total;
253+
254+
bsp_error_t error =
255+
matlab_to_bsp_array_allocator(mx_a, &out->values, bsp_matlab_allocator);
256+
if (error != BSP_SUCCESS) {
257+
mexErrMsgIdAndTxt("BinSparse:MemoryError",
258+
"Failed to allocate dense values");
259+
}
260+
}
261+
262+
void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
263+
if (nrhs < 1 || nrhs > 3) {
264+
mexErrMsgIdAndTxt(
265+
"BinSparse:InvalidArgs",
266+
"Usage: matrix = binsparse_from_ssmc(A [, Zeros] [, format])");
267+
}
268+
269+
if (nlhs > 1) {
270+
mexErrMsgIdAndTxt("BinSparse:TooManyOutputs", "Too many output arguments");
271+
}
272+
273+
const mxArray* mx_a = prhs[0];
274+
bsp_matrix_format_t target_format = parse_format(nrhs, prhs);
275+
276+
if (!mxIsSparse(mx_a)) {
277+
bsp_matrix_t dense_matrix;
278+
build_dense_matrix(mx_a, &dense_matrix, target_format);
279+
plhs[0] = bsp_matrix_to_matlab_struct(&dense_matrix);
280+
bsp_destroy_matrix_t(&dense_matrix);
281+
return;
282+
}
283+
284+
if (target_format == BSP_DMAT || target_format == BSP_DVEC) {
285+
mexErrMsgIdAndTxt("BinSparse:InvalidFormat",
286+
"Sparse matrix cannot use DMAT/DVEC formats");
287+
}
288+
289+
if (mxIsComplex(mx_a) || !mxIsDouble(mx_a)) {
290+
mexErrMsgIdAndTxt("BinSparse:InvalidMatrix",
291+
"A must be a real sparse double matrix");
292+
}
293+
294+
const mxArray* mx_zeros = NULL;
295+
if (nrhs >= 2 && mxIsSparse(prhs[1])) {
296+
mx_zeros = prhs[1];
297+
}
298+
299+
if (mx_zeros && (mxIsComplex(mx_zeros) || !mxIsDouble(mx_zeros))) {
300+
mexErrMsgIdAndTxt("BinSparse:InvalidZeros",
301+
"Zeros must be a real sparse double matrix");
302+
}
303+
304+
if (mx_zeros &&
305+
(mxGetM(mx_a) != mxGetM(mx_zeros) || mxGetN(mx_a) != mxGetN(mx_zeros))) {
306+
mexErrMsgIdAndTxt("BinSparse:DimensionMismatch",
307+
"A and Zeros must have matching dimensions");
308+
}
309+
310+
matlab_csc_t a_csc = {0};
311+
matlab_csc_t z_csc = {0};
312+
313+
if (extract_matlab_csc(mx_a, &a_csc) != 0) {
314+
mexErrMsgIdAndTxt("BinSparse:InvalidMatrix",
315+
"Failed to extract CSC data from A");
316+
}
317+
318+
bool have_zeros = false;
319+
if (mx_zeros) {
320+
if (extract_matlab_csc(mx_zeros, &z_csc) != 0) {
321+
mexErrMsgIdAndTxt("BinSparse:InvalidZeros",
322+
"Failed to extract CSC data from Zeros");
323+
}
324+
have_zeros = true;
325+
}
326+
327+
bsp_matrix_t csc_matrix;
328+
if (have_zeros) {
329+
build_csc_merged(&a_csc, &z_csc, &csc_matrix);
330+
} else {
331+
build_csc_from_a(&a_csc, &csc_matrix);
332+
}
333+
334+
bsp_matrix_t result = csc_matrix;
335+
if (target_format != BSP_CSC) {
336+
result = bsp_convert_matrix(csc_matrix, target_format);
337+
bsp_destroy_matrix_t(&csc_matrix);
338+
339+
if (result.format != target_format) {
340+
bsp_destroy_matrix_t(&result);
341+
mexErrMsgIdAndTxt("BinSparse:ConversionError",
342+
"Failed to convert matrix to requested format");
343+
}
344+
}
345+
346+
if (!matrix_uses_allocator(&result, bsp_matlab_allocator)) {
347+
bsp_matrix_t copied;
348+
bsp_error_t error =
349+
bsp_matrix_copy_with_allocator(&result, &copied, bsp_matlab_allocator);
350+
bsp_destroy_matrix_t(&result);
351+
if (error != BSP_SUCCESS) {
352+
mexErrMsgIdAndTxt("BinSparse:MemoryError",
353+
"Failed to allocate MATLAB-owned matrix");
354+
}
355+
result = copied;
356+
}
357+
358+
plhs[0] = bsp_matrix_to_matlab_struct(&result);
359+
bsp_destroy_matrix_t(&result);
360+
}

‎bindings/matlab/binsparse_from_ssmc.m‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
function matrix = binsparse_from_ssmc (A, Zeros, format)
2+
%BINSPARSE_FROM_SSMC convert SuiteSparse A+Zeros to a Binsparse matrix struct
3+
%
4+
% Usage:
5+
% matrix = binsparse_from_ssmc(A, Zeros)
6+
% matrix = binsparse_from_ssmc(A, Zeros, format)
7+
%
8+
% This is a thin wrapper over the binsparse_from_ssmc MEX function.
9+
10+
error('binsparse_from_ssmc mexFunction not found; compile with build_matlab_bindings first');

0 commit comments

Comments
 (0)