diff --git a/bindings/matlab/binsparse_to_ssmc_problem.m b/bindings/matlab/binsparse_to_ssmc_problem.m index 8c06695..0555ee9 100644 --- a/bindings/matlab/binsparse_to_ssmc_problem.m +++ b/bindings/matlab/binsparse_to_ssmc_problem.m @@ -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; @@ -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; @@ -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 @@ -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 @@ -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', ... @@ -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)); @@ -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]; @@ -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 @@ -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; @@ -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 diff --git a/bindings/matlab/test_binsparse_to_ssmc_problem.m b/bindings/matlab/test_binsparse_to_ssmc_problem.m index 3f998bf..0cbcddd 100644 --- a/bindings/matlab/test_binsparse_to_ssmc_problem.m +++ b/bindings/matlab/test_binsparse_to_ssmc_problem.m @@ -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);