Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
172 changes: 136 additions & 36 deletions bindings/matlab/binsparse_to_ssmc_problem.m
Original file line number Diff line number Diff line change
Expand Up @@ -143,12 +143,12 @@ function validate_problem(bsp_problem)

switch format
case 'DVEC'
values = matrix_values(bsp, count);
[values, iso] = matrix_values(bsp, count);
if n ~= 1 || count ~= m
error('BinSparse:InvalidMatrix', ...
'DVEC dimensions do not match the stored value count');
end
matrix = reshape(values, [m, 1]);
matrix = reshape(expand_iso(values, iso, count), [m, 1]);
Zeros = sparse(m, n);
return;

Expand All @@ -158,8 +158,8 @@ function validate_problem(bsp_problem)
error('BinSparse:InvalidMatrix', ...
'DMATR dimensions do not match the stored value count');
end
values = matrix_values(bsp, count);
matrix = reshape(values, [n, m]).';
[values, iso] = matrix_values(bsp, count);
matrix = reshape(expand_iso(values, iso, count), [n, m]).';
Zeros = sparse(m, n);
return;

Expand All @@ -169,24 +169,50 @@ function validate_problem(bsp_problem)
error('BinSparse:InvalidMatrix', ...
'DMATC dimensions do not match the stored value count');
end
values = matrix_values(bsp, count);
matrix = reshape(values, [m, n]);
[values, iso] = matrix_values(bsp, count);
matrix = reshape(expand_iso(values, iso, count), [m, n]);
Zeros = sparse(m, n);
return;
end

[rows, cols] = sparse_indices(bsp, format, m, n, count);
values = matrix_values(bsp, count);
[rows, cols, values] = expand_structure(bsp, rows, cols, values, m, n);
[values, iso] = matrix_values(bsp, count);
[rows, cols, values, iso] = ...
expand_structure(bsp, rows, cols, values, iso, m, n);

% An ISO matrix carries a single stored value, and sparse() expands a scalar
% over the index vectors, so its values are never materialized at nnz length.
% The entries are then either all zero or all nonzero, which also removes the
% need for a mask and for masked copies of the index vectors.
if iso
if split_zeros && values == 0
matrix = sparse(m, n);
Zeros = sparse(rows, cols, 1, m, n);
else
matrix = sparse(rows, cols, values, m, n);
Zeros = sparse(m, n);
end
return;
end

is_zero = (real(values) == 0) & (imag(values) == 0);
if split_zeros
matrix = sparse(rows(~is_zero), cols(~is_zero), ...
values(~is_zero), m, n);
Zeros = sparse(rows(is_zero), cols(is_zero), 1, m, n);
else
if ~split_zeros
matrix = sparse(rows, cols, values, m, n);
Zeros = sparse(m, n);
return;
end

% MATLAB's == compares both parts of a complex value, so this needs no
% separate imag() test; taking one would allocate an nnz-length array of
% zeros for the common case of real values.
is_zero = (values == 0);
if ~any(is_zero)
matrix = sparse(rows, cols, values, m, n);
Zeros = sparse(m, n);
else
keep = ~is_zero;
matrix = sparse(rows(keep), cols(keep), values(keep), m, n);
clear keep
Zeros = sparse(rows(is_zero), cols(is_zero), 1, m, n);
end
end

Expand Down Expand Up @@ -247,8 +273,13 @@ function validate_problem(bsp_problem)
end
end

function values = matrix_values(bsp, count)
function [values, iso] = matrix_values(bsp, count)
% Returns the stored values and whether they are still in ISO form, that is a
% single value standing for all count entries. The caller expands it only
% where an nnz-length array is genuinely required; the sparse constructors
% take the scalar directly.
values = bsp.values(:);
iso = false;
if logical(bsp.is_iso)
if count == 0
if numel(values) > 1
Expand All @@ -260,7 +291,7 @@ function validate_problem(bsp_problem)
error('BinSparse:InvalidMatrix', ...
'An ISO matrix must contain exactly one value');
else
values = repmat(values, count, 1);
iso = true;
end
elseif numel(values) ~= count
error('BinSparse:InvalidMatrix', ...
Expand All @@ -269,7 +300,14 @@ function validate_problem(bsp_problem)
values = exact_double(values);
end

function [rows, cols, values] = expand_structure(bsp, rows, cols, values, m, n)
function values = expand_iso(values, iso, count)
if iso
values = repmat(values, count, 1);
end
end

function [rows, cols, values, iso] = ...
expand_structure(bsp, rows, cols, values, iso, m, n)
structure = 'general';
if isfield(bsp, 'structure') && ~isempty(bsp.structure)
structure = lower(char(bsp.structure));
Expand All @@ -293,19 +331,33 @@ function validate_problem(bsp_problem)
'Stored entries do not match the declared matrix triangle');
end

off_diagonal = (rows ~= cols);
mirror_rows = cols(off_diagonal);
mirror_cols = rows(off_diagonal);
mirror_values = values(off_diagonal);
if starts_with(structure, 'hermitian_')
mirror_values = conj(mirror_values);
mirror = @conj;
elseif starts_with(structure, 'skew_symmetric_')
mirror_values = -mirror_values;
elseif ~starts_with(structure, 'symmetric_')
mirror = @uminus;
elseif starts_with(structure, 'symmetric_')
mirror = @(v) v;
else
error('BinSparse:UnsupportedStructure', ...
'Unsupported Binsparse matrix structure: %s', structure);
end

off_diagonal = (rows ~= cols);
mirror_rows = cols(off_diagonal);
mirror_cols = rows(off_diagonal);

% Mirroring an ISO matrix leaves it ISO whenever it leaves the single stored
% value alone, which covers every symmetric matrix and the real hermitian and
% all-zero skew-symmetric ones. Only the rest have to be expanded here.
if iso && isequaln(mirror(values), values)
rows = [rows; mirror_rows];
cols = [cols; mirror_cols];
return;
end

values = expand_iso(values, iso, numel(rows));
iso = false;
mirror_values = mirror(values(off_diagonal));
rows = [rows; mirror_rows];
cols = [cols; mirror_cols];
values = [values; mirror_values];
Expand Down Expand Up @@ -359,14 +411,32 @@ function validate_matrix_header(bsp)
if ~isnumeric(indices) || ~isreal(indices)
error('BinSparse:InvalidMatrix', '%s must be real numeric data', name);
end
indices = double(indices(:));
if ~isempty(expected) && numel(indices) ~= expected
error('BinSparse:InvalidMatrix', ...
'%s length does not match nnz', name);
end
if any(~isfinite(indices)) || any(indices < 0) || ...
any(fix(indices) ~= indices) || any(indices >= limit)
error('BinSparse:InvalidMatrix', '%s contains an invalid index', name);

if isinteger(indices)
% A Binsparse index array is normally stored as an integer type, where
% the finiteness and integrality tests are vacuous and the range tests
% reduce to min and max. Running them here, before the widening to
% double, keeps the check from building a second nnz-length temporary
% just to hold fix(indices).
if ~isempty(indices)
if double(max(indices(:))) >= limit || ...
(intmin(class(indices)) < 0 && double(min(indices(:))) < 0)
error('BinSparse:InvalidMatrix', ...
'%s contains an invalid index', name);
end
end
indices = double(indices(:));
else
indices = double(indices(:));
if any(~isfinite(indices)) || any(indices < 0) || ...
any(fix(indices) ~= indices) || any(indices >= limit)
error('BinSparse:InvalidMatrix', ...
'%s contains an invalid index', name);
end
end
end

Expand All @@ -393,17 +463,35 @@ function validate_matrix_header(bsp)
end

function require_ordered_pairs(first, second, format)
if numel(first) < 2
% Scanned in blocks, so the check costs a fixed few megabytes rather than
% several nnz-length temporaries. Consecutive blocks overlap by one entry so
% that no adjacent pair straddles a block boundary unchecked.
count = numel(first);
if count < 2
return;
end
bad = diff(first) < 0 | ...
(diff(first) == 0 & diff(second) <= 0);
if any(bad)
error('BinSparse:InvalidMatrix', ...
'%s indices are not sorted and unique', format);
lo = 1;
while lo < count
hi = min(count, lo + scan_block_size());
step = diff(first(lo:hi));
bad = step < 0;
tied = (step == 0);
clear step
if any(tied)
bad = bad | (tied & diff(second(lo:hi)) <= 0);
end
if any(bad)
error('BinSparse:InvalidMatrix', ...
'%s indices are not sorted and unique', format);
end
lo = hi;
end
end

function n = scan_block_size()
n = 4194304;
end

function require_segment_order(indices, pointers, format)
for k = 1:numel(pointers)-1
first = pointers(k) + 1;
Expand All @@ -423,8 +511,20 @@ function require_nonempty_segments(pointers, format)
end

function require_strictly_increasing(values, label)
if numel(values) > 1 && any(diff(values) <= 0)
error('BinSparse:InvalidMatrix', '%s are not sorted and unique', label);
% Blocked for the same reason as require_ordered_pairs: a CVEC index array is
% nnz long, so diff() over the whole of it is not affordable at scale.
count = numel(values);
if count < 2
return;
end
lo = 1;
while lo < count
hi = min(count, lo + scan_block_size());
if any(diff(values(lo:hi)) <= 0)
error('BinSparse:InvalidMatrix', ...
'%s are not sorted and unique', label);
end
lo = hi;
end
end

Expand Down
32 changes: 32 additions & 0 deletions bindings/matlab/test_binsparse_to_ssmc_problem.m
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,38 @@
assert(nnz(Problem.A) == 0);
assert(isequal(Problem.Zeros, sparse(rows + 1, cols + 1, 1, 3, 3)));

% An ISO matrix keeps its single stored value all the way to sparse(), which
% expands it, so these check that the expansion still lands on every entry and
% that mirroring drops the ISO form exactly when it changes the value.
raw = struct('metadata', metadata, ...
'A', make_matrix(4, rows, cols, [], 3, 3, 'COOR', true, 'general'));
Problem = binsparse_to_ssmc_problem(raw);
assert(isequal(Problem.A, sparse(rows + 1, cols + 1, 4, 3, 3)));
assert(~isfield(Problem, 'Zeros'), 'a nonzero ISO matrix has no Zeros');

iso_symmetric = make_matrix(7, [0; 1; 1], [0; 0; 1], ...
[], 2, 2, 'COO', true, 'symmetric_lower');
raw = struct('metadata', metadata, 'A', iso_symmetric);
Problem = binsparse_to_ssmc_problem(raw);
assert(isequal(Problem.A, sparse(7 * ones(2, 2))));

iso_skew = make_matrix(5, 1, 0, [], 2, 2, 'COO', true, ...
'skew_symmetric_lower');
raw = struct('metadata', metadata, 'A', iso_skew);
Problem = binsparse_to_ssmc_problem(raw);
assert(isequal(Problem.A, sparse([0 -5; 5 0])));

iso_hermitian = make_matrix(2 + 3i, 1, 0, [], 2, 2, 'COO', true, ...
'hermitian_lower');
raw = struct('metadata', metadata, 'A', iso_hermitian);
Problem = binsparse_to_ssmc_problem(raw);
assert(isequal(Problem.A, sparse([0 2-3i; 2+3i 0])));

raw = struct('metadata', metadata, 'A', formats{1}, ...
'b', make_matrix(9, [], [], [], 3, 1, 'DVEC', true, 'general'));
Problem = binsparse_to_ssmc_problem(raw);
assert(isequal(Problem.b, [9; 9; 9]), 'a dense ISO vector must expand');

lower = make_matrix([1; 2; 0; 3], [0; 1; 1; 2], [0; 0; 1; 2], ...
[], 3, 3, 'COO', false, 'symmetric_lower');
raw = struct('metadata', metadata, 'A', lower);
Expand Down
Loading