From 3203b92527dc89d4800bf965a142c45a761b0eef Mon Sep 17 00:00:00 2001 From: mmelnich Date: Thu, 8 Oct 2026 15:03:14 -0700 Subject: [PATCH] MATLAB bindings for CQRRPT and CQRRT Add randlapack.cqrrpt and randlapack.cqrrt next to randlapack.bqrrp, for tall matrices (m >= d_factor*n): - CQRRPT (rank-revealing QR with column pivoting): [Q, R, J] with the rank detected by CQRRPT as the inner dimension of Q and R; 'eps', 'd_factor' and 'state' options. Taken from the unmerged bindings-wave1 branch, with `rank` reset before the call (CQRRPT returns early on an all-zero input without setting it) and the J-zeroing comment updated (RandLAPACK now zeroes J itself; the MEX keeps doing it for older installs). - CQRRT (unpivoted Cholesky QR for full-rank input): [Q, R] like qr(A, 0), with Q factored in place into its output buffer; R = cqrrt(A) uses CQRRT's R-only mode and forms no Q. CQRRT does not detect rank: a rank-deficient input usually completes with a numerically singular R (documented and tested); an exactly singular sketch or a Cholesky failure is reported as randlapack:cqrrt_mex:breakdown, pointing at cqrrpt and bqrrp. - The RNG-state helpers move to +randlapack/private and bqrrp.m uses them. - matlab/CMakeLists.txt: randlapack_add_mex(name) applies the per-target setup (include dirs, output dir, static libstdc++ on GNU, Windows DLL staging and libiomp5md.dll removal, MSVC OpenMP linkage) to all three targets; bqrrp_mex's flags.make and link.txt are unchanged on Linux. - All three MEX files: errors raised by the MEX keep their identifier. feval("error") throws matlab::engine::MATLABExecutionException, which derives from matlab::engine::Exception and std::exception but not from matlab::Exception, so the old `catch (const matlab::Exception&)` never matched and the std::exception handler re-raised it as ...:StdError. A `catch (const matlab::engine::Exception&) { throw; }` now comes first. - Tests, demos and README for both drivers. The CQRRPT demo compares against qr(A, 'econ', 'vector') after a warm-up call (it used to compare against a full QR with an m-by-m Q). Help texts and README say "matches qr(A, 'econ', 'vector')" for bqrrp and cqrrpt. Verified on Linux (GCC 15.2, MATLAB R2026a, RandLAPACK main): test_bqrrp, test_cqrrpt and test_cqrrt pass, including rank detection, the all-zero and breakdown error IDs, R-only vs [Q, R] agreement, a rank-deficient input completing with a singular R, and a cond(A) = 1e10 input; CQRRT on 50000x500 matches qr(A, 0)'s R to 2e-15 (2.9x faster), CQRRPT on 20000x200 is 8.8x faster than qr(A, 'econ', 'vector'). On Windows (MSVC 19.44, MATLAB R2026a, RandLAPACK main installed with install.ps1, threaded oneMKL): the three MEX files build with the README command, import only MATLAB's libiomp5md for OpenMP, and all three test suites pass with KMP_DUPLICATE_LIB_OK unset. --- CMakeLists.txt | 2 +- README.md | 31 +- matlab/+randlapack/bqrrp.m | 36 +-- matlab/+randlapack/cqrrpt.m | 124 ++++++++ matlab/+randlapack/cqrrt.m | 109 +++++++ matlab/+randlapack/private/rl_default_state.m | 11 + .../+randlapack/private/rl_normalize_state.m | 28 ++ matlab/CMakeLists.txt | 149 +++++---- matlab/examples/cqrrpt_demo.m | 46 +++ matlab/examples/cqrrt_demo.m | 38 +++ matlab/src/bqrrp_mex.cc | 13 +- matlab/src/cqrrpt_mex.cc | 283 ++++++++++++++++++ matlab/src/cqrrt_mex.cc | 260 ++++++++++++++++ matlab/tests/test_cqrrpt.m | 181 +++++++++++ matlab/tests/test_cqrrt.m | 177 +++++++++++ 15 files changed, 1384 insertions(+), 104 deletions(-) create mode 100644 matlab/+randlapack/cqrrpt.m create mode 100644 matlab/+randlapack/cqrrt.m create mode 100644 matlab/+randlapack/private/rl_default_state.m create mode 100644 matlab/+randlapack/private/rl_normalize_state.m create mode 100644 matlab/examples/cqrrpt_demo.m create mode 100644 matlab/examples/cqrrt_demo.m create mode 100644 matlab/src/cqrrpt_mex.cc create mode 100644 matlab/src/cqrrt_mex.cc create mode 100644 matlab/tests/test_cqrrpt.m create mode 100644 matlab/tests/test_cqrrt.m diff --git a/CMakeLists.txt b/CMakeLists.txt index bdf6422..0dc8760 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -38,5 +38,5 @@ else() "Set -DMatlab_ROOT_DIR=/path/to/MATLAB to enable. " "The +randlapack package can still be used by adding " "${CMAKE_CURRENT_SOURCE_DIR}/matlab to the MATLAB path, " - "but calls into bqrrp_mex will fail until the MEX is built.") + "but calls into the MEX files will fail until they are built.") endif() diff --git a/README.md b/README.md index 8881bd3..2ce50b7 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,13 @@ MATLAB, and eventually Python, bindings for [RandLAPACK](https://github.com/Ball ## Scope (v0) -Proof-of-concept binding for a single driver, `BQRRP` (randomized blocked QR with column pivoting), in single and double precision, real matrices, column-major. More drivers as their C++ APIs stabilize. Python next. +MATLAB bindings for three QR drivers, in single and double precision, real matrices, column-major: + +* `BQRRP`: randomized blocked QR with column pivoting, any shape. +* `CQRRPT`: randomized rank-revealing QR with column pivoting, for tall matrices. +* `CQRRT`: randomized Cholesky QR without pivoting, for tall full-rank matrices. + +More drivers as their C++ APIs stabilize. Python next. ## Requirements @@ -24,7 +30,7 @@ cmake --build build -j Pass `-DCMAKE_BUILD_TYPE=Release`: without a build type the MEX compiles without optimization (Ninja with MSVC even defaults to Debug). -If MATLAB is not found, MEX targets are skipped with a clear warning; the `.m` files still ship and can be used once the MEX is built. On success, `matlab/+randlapack/private/bqrrp_mex.` is in place. +If MATLAB is not found, MEX targets are skipped with a clear warning; the `.m` files still ship and can be used once the MEX is built. On success, `bqrrp_mex.`, `cqrrpt_mex.` and `cqrrt_mex.` are in `matlab/+randlapack/private/`. ### Windows @@ -56,12 +62,29 @@ err = norm(A(:, J) - Q*R, 'fro') / norm(A, 'fro'); Two output modes, chosen by a trailing string (qr-style): ```matlab -[Q, R, J] = randlapack.bqrrp(A, 'explicit') % default; matches qr(A, 'vector') +[Q, R, J] = randlapack.bqrrp(A, 'explicit') % default; matches qr(A, 'econ', 'vector') [A_out, tau, J] = randlapack.bqrrp(A, 'implicit') % GEQP3-format; skips Q materialization ``` Run `help randlapack.bqrrp` for the full signature, including the optional `b_sz`, `d_factor`, and RNG `state` arguments. +For **tall** matrices (`m >= d_factor*n`, with `d_factor = 1.25` by default), two sketch-and-precondition drivers: + +```matlab +A = randn(20000, 200); + +% CQRRPT: rank-revealing, pivoted. The detected rank is size(R, 1). +[Q, R, J] = randlapack.cqrrpt(A); +err = norm(A(:, J) - Q*R, 'fro') / norm(A, 'fro'); + +% CQRRT: unpivoted, for full-rank A. Same layout as qr(A, 0). +[Q, R] = randlapack.cqrrt(A); +err = norm(A - Q*R, 'fro') / norm(A, 'fro'); +R = randlapack.cqrrt(A); % R only: Q is not formed +``` + +CQRRT does not detect rank: on a rank-deficient input it usually returns without error, with a numerically singular R (it stops with an error only on an exactly singular sketch, e.g. a zero column, or a Cholesky failure). If the rank is uncertain, check `abs(diag(R))` or use `randlapack.cqrrpt`. Both take `'d_factor'` and `'state'` name-value options (CQRRPT also `'eps'`); see `help randlapack.cqrrpt` and `help randlapack.cqrrt`. + ## Runtime notes On **Linux MATLAB R2026a + Ubuntu 24.04** (and similar new-glibc hosts), MATLAB ships an older `libstdc++.so.6` than your compiler's. MEX files built against the system `libstdc++` fail at load with `version GLIBCXX_3.4.32 not found`. Workaround: @@ -76,5 +99,5 @@ On **Windows**, the MEX uses MATLAB's own OpenMP runtime (`libiomp5md`) instead ## Conventions -* `J` is **1-based**, matching MATLAB. BQRRP stores 1-based pivots natively; no conversion in the MEX layer. +* `J` is **1-based**, matching MATLAB. BQRRP and CQRRPT store 1-based pivots natively; no conversion in the MEX layer. * RNG `state` is Philox4x32: `state.counter` is `uint32[4]`, `state.key` is `uint32[2]`. Omit it (defaults to seed 0) or pass a scalar `uint32` seed. diff --git a/matlab/+randlapack/bqrrp.m b/matlab/+randlapack/bqrrp.m index 2684331..56a21b0 100644 --- a/matlab/+randlapack/bqrrp.m +++ b/matlab/+randlapack/bqrrp.m @@ -3,7 +3,7 @@ % % Two output modes, selected via a trailing string argument: % -% 'explicit' (default): matches MATLAB qr(A, 'vector') +% 'explicit' (default): matches MATLAB qr(A, 'econ', 'vector') % [Q, R, J] = randlapack.bqrrp(A) % [Q, R, J] = randlapack.bqrrp(A, 'explicit') % Q is m-by-k orthogonal, R is k-by-n upper-triangular, k = min(m, n), @@ -50,7 +50,7 @@ % pass a smaller value to actually exercise the blocked structure. b_sz = int64(min(m, n)); d_factor = cast(1.25, class(A)); - state = make_default_state(); + state = rl_default_state(); mode = 'explicit'; % --- Argument parsing --- @@ -75,7 +75,7 @@ d_factor = cast(positional{2}, class(A)); end if numel(positional) >= 3 && ~isempty(positional{3}) - state = normalize_state(positional{3}); + state = rl_normalize_state(positional{3}, mfilename); end if numel(positional) >= 4 error('randlapack:bqrrp:nargin', ... @@ -104,33 +104,3 @@ % out2 = tau (length-n Householder scalars) [out1, out2, J, state_out] = bqrrp_mex(A, b_sz, d_factor, state, mode); end - - -function state = make_default_state() - state.counter = uint32([0, 0, 0, 0]); - state.key = uint32([0, 0]); -end - - -function state = normalize_state(in) - if isstruct(in) - if ~isfield(in, 'counter') || ~isfield(in, 'key') - error('randlapack:bqrrp:state', ... - 'state struct must have fields .counter and .key'); - end - validateattributes(in.counter, {'uint32'}, {'vector', 'numel', 4}, ... - mfilename, 'state.counter'); - validateattributes(in.key, {'uint32'}, {'vector', 'numel', 2}, ... - mfilename, 'state.key'); - state.counter = uint32(in.counter(:)).'; - state.key = uint32(in.key(:)).'; - elseif isnumeric(in) && isscalar(in) - % Treat scalar as a seed; place in key[1]. Counter starts at zero. - seed = uint32(in); - state.counter = uint32([0, 0, 0, 0]); - state.key = uint32([0, seed]); - else - error('randlapack:bqrrp:state', ... - 'state must be a struct or a scalar seed'); - end -end diff --git a/matlab/+randlapack/cqrrpt.m b/matlab/+randlapack/cqrrpt.m new file mode 100644 index 0000000..03a5318 --- /dev/null +++ b/matlab/+randlapack/cqrrpt.m @@ -0,0 +1,124 @@ +function [Q, R, J, state_out] = cqrrpt(A, varargin) +%RANDLAPACK.CQRRPT Randomized rank-revealing QR with column pivoting. +% +% [Q, R, J] = randlapack.cqrrpt(A) computes a rank-revealing QR +% factorization of a TALL matrix A, so that +% A(:, J) ~= Q * R +% with Q m-by-k orthonormal, R k-by-n upper-trapezoidal, and J a length-n +% vector of 1-based column pivots. Output layout matches +% qr(A, 'econ', 'vector'). +% +% k is DISCOVERED, not requested: CQRRPT reveals the numerical rank of A +% and that rank is the inner dimension of Q and R. On a full-rank input +% k = n. Read it off size(R, 1). +% +% What the factorization guarantees, and what it does not: Q is +% orthonormal and the LEADING k pivoted columns always reproduce exactly, +% A(:, J(1:k)) == Q * R(:, 1:k) +% to machine precision. The full identity A(:, J) == Q*R additionally +% needs k to have reached the numerical rank of A: the ordinary case, +% including on rank-deficient input. It can fall short when CQRRPT's +% internal Cholesky QR breaks down and the conservative fallback estimate +% takes over, so check size(R, 1) and the residual on inputs whose +% conditioning you do not know. +% +% CQRRPT sketches A down to d = d_factor*n rows, so it requires +% m >= d_factor*n. For wide or square matrices use randlapack.bqrrp. +% +% Name-value options: +% 'eps' conditioning tolerance for rank detection. Default: +% sqrt(eps(class(A))). +% +% eps is NOT a singular-value truncation threshold. +% CQRRPT first truncates using a fixed machine-epsilon +% cutoff, so columns well above machine precision are +% kept whatever eps says. eps takes effect only in the +% re-estimation that runs when the internal Cholesky QR +% fails (e.g. on an ill-conditioned input), where the rank +% is cut once the R diagonal spans a ratio of +% sqrt(eps/eps(class(A))). It therefore bounds the +% orthogonality loss, which for Cholesky QR scales like +% u*cond(R)^2. Smaller eps means a stricter bound and a +% smaller returned rank. +% 'd_factor' sketch embedding factor (>= 1.0). Default: 1.25. +% Larger values sketch more rows: steadier rank +% detection, more work, and a taller input requirement. +% 'state' RNG state: a uint32 scalar seed, or a struct with +% .counter (uint32[4]) and .key (uint32[2]). +% +% [..., state] = randlapack.cqrrpt(...) returns the advanced RNG state. +% +% A must be single or double precision, real, and 2D. +% +% Example: +% B = randn(5000, 40); +% A = [B, B(:, 1:10)]; % 50 columns, rank 40 +% [Q, R, J] = randlapack.cqrrpt(A); +% detected_rank = size(R, 1) % ~40 +% err = norm(A(:, J) - Q*R, 'fro') / norm(A, 'fro'); +% +% See also QR, RANDLAPACK.BQRRP, RANDLAPACK.CQRRT. + + % --- Required argument check --- + validateattributes(A, {'single', 'double'}, {'2d', 'real', 'finite'}, ... + mfilename, 'A', 1); + [m, n] = size(A); + if m == 0 || n == 0 + error('randlapack:cqrrpt:emptyInput', 'A must be nonempty'); + end + + % --- Defaults --- + % eps default: sqrt of machine epsilon in the working precision. eps + % bounds the orthogonality loss in the fallback rank estimate (see the + % help above), and Cholesky QR loses orthogonality like u*cond(R)^2, so + % sqrt(u) is its natural scale. + opts = struct('eps', sqrt(eps(class(A))), ... + 'd_factor', 1.25, ... + 'state', rl_default_state()); + + % --- Name-value parsing --- + if mod(numel(varargin), 2) ~= 0 + error('randlapack:cqrrpt:nargin', ... + 'Options must be given as name-value pairs.'); + end + known = fieldnames(opts); + for i = 1:2:numel(varargin) + name = varargin{i}; + if ~(ischar(name) || (isstring(name) && isscalar(name))) + error('randlapack:cqrrpt:optionName', ... + 'Option names must be character vectors or strings.'); + end + name = char(name); + idx = find(strcmpi(known, name), 1); + if isempty(idx) + error('randlapack:cqrrpt:unknownOption', ... + 'Unknown option ''%s''. Valid options: %s.', ... + name, strjoin(known', ', ')); + end + opts.(known{idx}) = varargin{i + 1}; + end + + % --- Validation --- + validateattributes(opts.eps, {'numeric'}, {'scalar', 'real', 'nonnegative'}, ... + mfilename, 'eps'); + validateattributes(opts.d_factor, {'numeric'}, {'scalar', 'real', '>=', 1}, ... + mfilename, 'd_factor'); + + % The tall-matrix requirement is checked here as well as in the MEX, so + % the common mistake is caught before any marshalling happens and the + % message can name the MATLAB-level alternative. + d = floor(opts.d_factor * n); + if m < d + error('randlapack:cqrrpt:notTall', ... + ['CQRRPT requires m >= d_factor*n. Got m=%d, n=%d, ' ... + 'd_factor=%g (needs m >= %d). Lower d_factor toward 1.0, ' ... + 'or use randlapack.bqrrp, which has no tall requirement.'], ... + m, n, opts.d_factor, d); + end + + state = rl_normalize_state(opts.state, mfilename); + + % --- Dispatch to MEX --- + [Q, R, J, state_out] = cqrrpt_mex(A, cast(opts.d_factor, class(A)), ... + cast(opts.eps, class(A)), state); +end diff --git a/matlab/+randlapack/cqrrt.m b/matlab/+randlapack/cqrrt.m new file mode 100644 index 0000000..1ac1f88 --- /dev/null +++ b/matlab/+randlapack/cqrrt.m @@ -0,0 +1,109 @@ +function [out1, R, state_out] = cqrrt(A, varargin) +%RANDLAPACK.CQRRT Randomized Cholesky QR of a tall full-rank matrix. +% +% [Q, R] = randlapack.cqrrt(A) computes an unpivoted QR factorization of a +% TALL, FULL-RANK matrix A, so that +% A ~= Q * R +% with Q m-by-n orthonormal and R n-by-n upper-triangular. Output layout +% matches qr(A, 0) (the economy-size QR). +% +% R = randlapack.cqrrt(A), with a single output, returns R only and skips +% forming Q (CQRRT's R-only mode), like R = qr(A, 0). This is the cheaper +% call when only R is needed, e.g. to solve least-squares problems or to +% precondition with R. +% +% How it works: CQRRT sketches A down to d = d_factor*n rows, takes a QR of +% the small sketch, uses its R factor to precondition A, and finishes with +% Cholesky QR on the well-conditioned result. It does not pivot and does +% NOT detect rank: on a rank-deficient A it usually returns without error, +% with Q*R still reproducing A but R numerically singular (tiny diagonal +% entries). It stops with an error only when the sketch's R is exactly +% singular (for example, A has a zero column) or the Cholesky step fails. +% If the rank of A is uncertain, check abs(diag(R)) or use +% randlapack.cqrrpt, which reveals the rank (or randlapack.bqrrp). +% +% CQRRT requires m >= d_factor*n. For wide or square matrices use +% randlapack.bqrrp. +% +% Name-value options: +% 'd_factor' sketch embedding factor (>= 1.0). Default: 1.25. +% Larger values sketch more rows: a better +% preconditioner, more work, and a taller input +% requirement. +% 'state' RNG state: a uint32 scalar seed, or a struct with +% .counter (uint32[4]) and .key (uint32[2]). +% +% [Q, R, state] = randlapack.cqrrt(...) returns the advanced RNG state. +% (The single-output R-only form does not return it.) +% +% A must be single or double precision, real, and 2D. +% +% Example: +% A = randn(20000, 200); +% [Q, R] = randlapack.cqrrt(A); +% err = norm(A - Q*R, 'fro') / norm(A, 'fro'); +% R_only = randlapack.cqrrt(A); % no Q formed +% +% See also QR, RANDLAPACK.CQRRPT, RANDLAPACK.BQRRP. + + % --- Required argument check --- + validateattributes(A, {'single', 'double'}, {'2d', 'real', 'finite'}, ... + mfilename, 'A', 1); + [m, n] = size(A); + if m == 0 || n == 0 + error('randlapack:cqrrt:emptyInput', 'A must be nonempty'); + end + + % --- Defaults --- + opts = struct('d_factor', 1.25, ... + 'state', rl_default_state()); + + % --- Name-value parsing --- + if mod(numel(varargin), 2) ~= 0 + error('randlapack:cqrrt:nargin', ... + 'Options must be given as name-value pairs.'); + end + known = fieldnames(opts); + for i = 1:2:numel(varargin) + name = varargin{i}; + if ~(ischar(name) || (isstring(name) && isscalar(name))) + error('randlapack:cqrrt:optionName', ... + 'Option names must be character vectors or strings.'); + end + name = char(name); + idx = find(strcmpi(known, name), 1); + if isempty(idx) + error('randlapack:cqrrt:unknownOption', ... + 'Unknown option ''%s''. Valid options: %s.', ... + name, strjoin(known', ', ')); + end + opts.(known{idx}) = varargin{i + 1}; + end + + % --- Validation --- + validateattributes(opts.d_factor, {'numeric'}, {'scalar', 'real', '>=', 1}, ... + mfilename, 'd_factor'); + + % The tall-matrix requirement is checked here as well as in the MEX, so + % the common mistake is caught before any marshalling happens and the + % message can name the MATLAB-level alternative. + d = floor(opts.d_factor * n); + if m < d + error('randlapack:cqrrt:notTall', ... + ['CQRRT requires m >= d_factor*n. Got m=%d, n=%d, ' ... + 'd_factor=%g (needs m >= %d). Lower d_factor toward 1.0, ' ... + 'or use randlapack.bqrrp, which has no tall requirement.'], ... + m, n, opts.d_factor, d); + end + + state = rl_normalize_state(opts.state, mfilename); + + % --- Dispatch to MEX --- + d_factor = cast(opts.d_factor, class(A)); + if nargout <= 1 + % R-only mode: Q is not formed. + [~, out1] = cqrrt_mex(A, d_factor, false, state); + else + [out1, R, state_out] = cqrrt_mex(A, d_factor, true, state); + end +end diff --git a/matlab/+randlapack/private/rl_default_state.m b/matlab/+randlapack/private/rl_default_state.m new file mode 100644 index 0000000..710f298 --- /dev/null +++ b/matlab/+randlapack/private/rl_default_state.m @@ -0,0 +1,11 @@ +function state = rl_default_state() +%RL_DEFAULT_STATE Default Philox4x32 RNG state (all zeros). +% +% Shared by the randlapack.* wrappers. Functions in a private/ directory +% are visible to the enclosing package's functions and to nobody else, +% which is exactly the scope these helpers want. The compiled MEX files +% land in this same directory, for the same reason. + + state.counter = uint32([0, 0, 0, 0]); + state.key = uint32([0, 0]); +end diff --git a/matlab/+randlapack/private/rl_normalize_state.m b/matlab/+randlapack/private/rl_normalize_state.m new file mode 100644 index 0000000..e67fafa --- /dev/null +++ b/matlab/+randlapack/private/rl_normalize_state.m @@ -0,0 +1,28 @@ +function state = rl_normalize_state(in, caller) +%RL_NORMALIZE_STATE Coerce a user-supplied RNG state into the MEX form. +% +% Accepts a struct with .counter (uint32[4]) and .key (uint32[2]), or a +% scalar numeric seed (placed in key(2), with the counter left at zero). +% CALLER names the user-facing function, so the raised error identifier +% and message point at what the user actually called. + + if isstruct(in) + if ~isfield(in, 'counter') || ~isfield(in, 'key') + error(sprintf('randlapack:%s:state', caller), ... + 'state struct must have fields .counter and .key'); + end + validateattributes(in.counter, {'uint32'}, {'vector', 'numel', 4}, ... + caller, 'state.counter'); + validateattributes(in.key, {'uint32'}, {'vector', 'numel', 2}, ... + caller, 'state.key'); + state.counter = uint32(in.counter(:)).'; + state.key = uint32(in.key(:)).'; + elseif isnumeric(in) && isscalar(in) + % Treat a scalar as a seed; place it in key(2). Counter starts at zero. + state.counter = uint32([0, 0, 0, 0]); + state.key = uint32([0, uint32(in)]); + else + error(sprintf('randlapack:%s:state', caller), ... + 'state must be a struct or a scalar seed'); + end +end diff --git a/matlab/CMakeLists.txt b/matlab/CMakeLists.txt index 3de2331..93ec584 100644 --- a/matlab/CMakeLists.txt +++ b/matlab/CMakeLists.txt @@ -12,7 +12,7 @@ file(MAKE_DIRECTORY ${MEX_OUTPUT_DIR}) # symbols required by code built with GCC 13+. Statically linking avoids # the runtime dependency entirely -- the MEX file ships its own libstdc++ # symbols and is loadable in any MATLAB without LD_PRELOAD gymnastics. -# Applied via target_link_options on each matlab_add_mex target below. +# Applied via target_link_options in randlapack_add_mex below. # # These are GCC/Linux-only flags. AppleClang on macOS rejects -static-libgcc # ("unsupported option"): there is no libgcc to statically link (macOS uses @@ -24,68 +24,89 @@ else() set(MEX_STATIC_LIBSTDCXX_FLAGS "") endif() -matlab_add_mex( - NAME bqrrp_mex - SRC src/bqrrp_mex.cc - LINK_TO RandLAPACK::RandLAPACK - R2018a -) -# matlab_add_mex's LINK_TO doesn't propagate INTERFACE_INCLUDE_DIRECTORIES from -# the linked targets in all CMake versions. Expose them explicitly. We use -# target_include_directories (not a second target_link_libraries) because -# FindMatlab itself uses the plain (no-keyword) target_link_libraries form -# internally and CMake disallows mixing keyword/plain calls on one target. -target_include_directories(bqrrp_mex PRIVATE - $ -) -set_target_properties(bqrrp_mex PROPERTIES - LIBRARY_OUTPUT_DIRECTORY ${MEX_OUTPUT_DIR} - RUNTIME_OUTPUT_DIRECTORY ${MEX_OUTPUT_DIR} -) -target_compile_features(bqrrp_mex PRIVATE cxx_std_20) -target_link_options(bqrrp_mex PRIVATE ${MEX_STATIC_LIBSTDCXX_FLAGS}) -# Windows: stage the BLAS backend's runtime DLLs next to the MEX. MATLAB loads -# a MEX with the MEX file's own directory on the DLL search path, so app-local -# staging (the idiomatic Windows layout RandLAPACK's installer uses for its own -# executables) makes the MEX loadable without PATH edits. The helper is a -# no-op off Windows and only exists in RandLAPACK installs that export it. -if(COMMAND randlapack_stage_runtime_dlls) - randlapack_stage_runtime_dlls(bqrrp_mex) - # Never leave Intel's OpenMP runtime, libiomp5md.dll, beside the MEX. - # MATLAB already hosts it, and a second copy loaded from this directory - # aborts with OMP Error #15 at the first parallel region; this happens - # when the BLAS backend is threaded through Intel OpenMP and the RandLAPACK - # install records that runtime's directory among the DLL directories to - # stage. Without the copy, the backend binds to MATLAB's runtime. CMake - # runs a target's POST_BUILD commands in the order they are added (every - # generator does; the docs promise order only within one call), so this - # follows the staging. - if(WIN32) - add_custom_command(TARGET bqrrp_mex POST_BUILD - COMMAND ${CMAKE_COMMAND} -E rm -f - $/libiomp5md.dll - VERBATIM) +# Defines one MEX target, built from src/.cc, with everything every +# target here needs. Each driver gets one call at the bottom of this file. +function(randlapack_add_mex name) + matlab_add_mex( + NAME ${name} + SRC src/${name}.cc + LINK_TO RandLAPACK::RandLAPACK + R2018a + ) + # matlab_add_mex's LINK_TO doesn't propagate INTERFACE_INCLUDE_DIRECTORIES + # from the linked targets in all CMake versions. Expose them explicitly. We + # use target_include_directories (not a second target_link_libraries) + # because FindMatlab itself uses the plain (no-keyword) + # target_link_libraries form internally and CMake disallows mixing + # keyword/plain calls on one target. + target_include_directories(${name} PRIVATE + $ + ) + set_target_properties(${name} PROPERTIES + LIBRARY_OUTPUT_DIRECTORY ${MEX_OUTPUT_DIR} + RUNTIME_OUTPUT_DIRECTORY ${MEX_OUTPUT_DIR} + ) + target_compile_features(${name} PRIVATE cxx_std_20) + target_link_options(${name} PRIVATE ${MEX_STATIC_LIBSTDCXX_FLAGS}) + + # Windows: stage the BLAS backend's runtime DLLs next to the MEX. MATLAB + # loads a MEX with the MEX file's own directory on the DLL search path, so + # app-local staging (the idiomatic Windows layout RandLAPACK's installer + # uses for its own executables) makes the MEX loadable without PATH edits. + # The helper is a no-op off Windows and only exists in RandLAPACK installs + # that export it. + if(COMMAND randlapack_stage_runtime_dlls) + randlapack_stage_runtime_dlls(${name}) + # Never leave Intel's OpenMP runtime, libiomp5md.dll, beside the MEX. + # MATLAB already hosts it, and a second copy loaded from this directory + # aborts with OMP Error #15 at the first parallel region; this happens + # when the BLAS backend is threaded through Intel OpenMP and the + # RandLAPACK install records that runtime's directory among the DLL + # directories to stage. Without the copy, the backend binds to MATLAB's + # runtime. CMake runs a target's POST_BUILD commands in the order they + # are added (every generator does; the docs promise order only within + # one call), so this follows the staging. + if(WIN32) + add_custom_command(TARGET ${name} POST_BUILD + COMMAND ${CMAKE_COMMAND} -E rm -f + $/libiomp5md.dll + VERBATIM) + # Every MEX target stages the same DLL set (about 1 GB for + # oneMKL) into the same directory. Run in parallel, two targets' + # copies, or one target's copy and another's removal above, can + # reach the same file at once and fail with a sharing violation, + # so each MEX target waits for the one defined before it. + get_property(previous_mex GLOBAL PROPERTY RANDLAPACK_BINDINGS_LAST_MEX) + if(previous_mex) + add_dependencies(${name} ${previous_mex}) + endif() + set_property(GLOBAL PROPERTY RANDLAPACK_BINDINGS_LAST_MEX ${name}) + endif() endif() -endif() -# Windows (MSVC): MATLAB's process already hosts Intel's OpenMP runtime -# (libiomp5md), while code compiled with /openmp:llvm requests MSVC's runtime -# (default library libomp.lib, which loads libomp140.x86_64.dll); two OpenMP -# runtimes in one process abort at the first parallel region (OMP Error #15). -# Intel's runtime and LLVM's share the kmpc entry points (LLVM's libomp -# descends from Intel's), so resolve the OpenMP references from MATLAB's own -# import library (shipped in bin/win64 for exactly this purpose) and drop the -# compiler's default runtime (libompd.lib is its /MDd variant). MATLAB's -# runtime lacks some newer entry points (e.g. the ones MSVC emits for -# collapse, used in rl_rpchol.hh); code that needs them then fails to link, -# or, when a threaded BLAS++ brings Intel's newer import library, fails to -# load ("The specified procedure could not be found"), instead of silently -# loading a second runtime. Plain target_link_libraries form, matching -# FindMatlab's internal usage (keyword/plain forms cannot mix). Guarded on -# MSVC, not WIN32: these are MSVC linker flags, and MinGW builds use their own -# OpenMP runtime. -if(MSVC) - target_link_libraries(bqrrp_mex "${Matlab_ROOT_DIR}/bin/win64/libiomp5md.lib") - target_link_options(bqrrp_mex PRIVATE - "/nodefaultlib:libomp.lib" "/nodefaultlib:libompd.lib") -endif() + # Windows (MSVC): MATLAB's process already hosts Intel's OpenMP runtime + # (libiomp5md), while code compiled with /openmp:llvm requests MSVC's + # runtime (default library libomp.lib, which loads libomp140.x86_64.dll); + # two OpenMP runtimes in one process abort at the first parallel region + # (OMP Error #15). Intel's runtime and LLVM's share the kmpc entry points + # (LLVM's libomp descends from Intel's), so resolve the OpenMP references + # from MATLAB's own import library (shipped in bin/win64 for exactly this + # purpose) and drop the compiler's default runtime (libompd.lib is its + # /MDd variant). MATLAB's runtime lacks some newer entry points (e.g. the + # ones MSVC emits for collapse, used in rl_rpchol.hh); code that needs them + # then fails to link, or, when a threaded BLAS++ brings Intel's newer + # import library, fails to load ("The specified procedure could not be + # found"), instead of silently loading a second runtime. Plain + # target_link_libraries form, matching FindMatlab's internal usage + # (keyword/plain forms cannot mix). Guarded on MSVC, not WIN32: these are + # MSVC linker flags, and MinGW builds use their own OpenMP runtime. + if(MSVC) + target_link_libraries(${name} "${Matlab_ROOT_DIR}/bin/win64/libiomp5md.lib") + target_link_options(${name} PRIVATE + "/nodefaultlib:libomp.lib" "/nodefaultlib:libompd.lib") + endif() +endfunction() + +randlapack_add_mex(bqrrp_mex) +randlapack_add_mex(cqrrpt_mex) +randlapack_add_mex(cqrrt_mex) diff --git a/matlab/examples/cqrrpt_demo.m b/matlab/examples/cqrrpt_demo.m new file mode 100644 index 0000000..c542972 --- /dev/null +++ b/matlab/examples/cqrrpt_demo.m @@ -0,0 +1,46 @@ +% cqrrpt_demo Rank revelation with RandLAPACK CQRRPT from MATLAB. +% +% CQRRPT discovers the numerical rank of a tall matrix rather than taking +% it as a parameter. This demo builds a matrix of known rank by repeating +% columns, then checks that the detected rank, the pivots, and the residual +% all agree with the construction. +% +% The timings compare against MATLAB's economy-size column-pivoted QR, +% qr(A, 'econ', 'vector'), after a warm-up call that loads the MEX file. The +% randomized method pays off as the matrix grows. +% +% Run from MATLAB after: +% addpath('/path/to/randlapack-bindings/matlab') + +m = 20000; r = 160; extra = 40; +B = randn(m, r); +A = [B, B(:, 1:extra)]; % rank r, r+extra columns +perm = randperm(size(A, 2)); +A = A(:, perm); % scatter the duplicates + +randlapack.cqrrpt(randn(500, 50)); % warm-up: loads the MEX file + +t = tic; +[Q, R, J] = randlapack.cqrrpt(A); +dt = toc(t); + +k = size(R, 1); +fprintf('A is %d-by-%d, constructed rank %d\n', m, size(A, 2), r); +fprintf('CQRRPT detected rank %d in %.3f s\n', k, dt); +fprintf('||A(:, J) - Q*R||_F / ||A||_F = %.2e\n', ... + norm(A(:, J) - Q * R, 'fro') / norm(A, 'fro')); +fprintf('||Q^T Q - I||_F = %.2e\n', ... + norm(Q' * Q - eye(k), 'fro')); + +% The first k pivots should select an independent set. Checking against +% MATLAB's own rank() is the honest comparison. +fprintf('rank(A(:, J(1:k))) = %d (want %d)\n', ... + rank(A(:, J(1:k))), k); + +% Compare against MATLAB's dense QRCP on the same matrix (economy size, so +% that, like CQRRPT, it forms no m-by-m Q). +t = tic; [~, ~, Jd] = qr(A, 'econ', 'vector'); dt_dense = toc(t); +fprintf('\ndense qr(A, ''econ'', ''vector''): %.3f s\n', dt_dense); +fprintf('CQRRPT speedup: %.1fx\n', dt_dense / dt); +fprintf('leading pivots agree on %d of the first %d columns\n', ... + numel(intersect(J(1:k), Jd(1:k))), k); diff --git a/matlab/examples/cqrrt_demo.m b/matlab/examples/cqrrt_demo.m new file mode 100644 index 0000000..39bf9c4 --- /dev/null +++ b/matlab/examples/cqrrt_demo.m @@ -0,0 +1,38 @@ +% cqrrt_demo Tall unpivoted QR with RandLAPACK CQRRT from MATLAB. +% +% CQRRT factors a tall, full-rank matrix A = Q*R by sketching A, using the +% sketch's R to precondition A, and finishing with Cholesky QR. This demo +% compares it with MATLAB's economy-size qr(A, 0), then shows the R-only +% call, which skips forming Q. +% +% CQRRT expects full rank. For rank-deficient or very ill-conditioned +% inputs use randlapack.cqrrpt, which reveals the rank. +% +% Run from MATLAB after: +% addpath('/path/to/randlapack-bindings/matlab') + +m = 50000; n = 500; +A = randn(m, n); +randlapack.cqrrt(randn(500, 50)); % warm-up: loads the MEX file + +t = tic; +[Q, R] = randlapack.cqrrt(A); +dt = toc(t); +fprintf('A is %d-by-%d\n', m, n); +fprintf('CQRRT [Q, R]: %.3f s\n', dt); +fprintf('||A - Q*R||_F / ||A||_F = %.2e\n', norm(A - Q * R, 'fro') / norm(A, 'fro')); +fprintf('||Q^T Q - I||_F = %.2e\n', norm(Q' * Q - eye(n), 'fro')); + +t = tic; [Qd, Rd] = qr(A, 0); dt_dense = toc(t); +fprintf('\ndense qr(A, 0): %.3f s\n', dt_dense); +fprintf('CQRRT speedup: %.1fx\n', dt_dense / dt); + +% R is unique up to the signs of its rows, so compare |R| entries. +fprintf('max | |R| - |R_qr| | / max|R| = %.2e\n', ... + max(abs(abs(R(:)) - abs(Rd(:)))) / max(abs(R(:)))); + +% R-only: the cheaper call when Q is not needed (least squares, preconditioning). +t = tic; +R_only = randlapack.cqrrt(A); +dt_r = toc(t); +fprintf('\nCQRRT R only: %.3f s\n', dt_r); diff --git a/matlab/src/bqrrp_mex.cc b/matlab/src/bqrrp_mex.cc index aee4951..f0787e2 100644 --- a/matlab/src/bqrrp_mex.cc +++ b/matlab/src/bqrrp_mex.cc @@ -284,9 +284,18 @@ class MexFunction : public matlab::mex::Function { raise("randlapack:bqrrp_mex:RandLAPACKError", e.what()); } catch (const RandBLAS::Error& e) { raise("randlapack:bqrrp_mex:RandBLASError", e.what()); + } catch (const matlab::engine::Exception&) { + // Errors raised by this MEX (raise() calls feval("error"), which + // throws matlab::engine::MATLABException) and other MATLAB engine + // errors propagate unchanged, so MATLAB sees the original error ID + // and message. In the MEX API (cppmex/mexException.hpp) these + // derive from std::exception but not from matlab::Exception (the + // C++ Engine API's classes of the same name differ): without this + // handler, the std::exception one below would re-raise them as + // ...:StdError. + throw; } catch (const matlab::Exception&) { - // MATLAB exceptions (typically from feval(error)) propagate through - // unchanged so MATLAB sees the original error ID and message. + // MATLAB Data API exceptions propagate unchanged. throw; } catch (const std::exception& e) { raise("randlapack:bqrrp_mex:StdError", e.what()); diff --git a/matlab/src/cqrrpt_mex.cc b/matlab/src/cqrrpt_mex.cc new file mode 100644 index 0000000..85081e6 --- /dev/null +++ b/matlab/src/cqrrpt_mex.cc @@ -0,0 +1,283 @@ +// MEX wrapper for RandLAPACK::CQRRPT using MATLAB's C++ Data API (R2018a+). +// +// MATLAB-side signature (low-level; users normally call randlapack.cqrrpt.m): +// +// [Q, R, J, state_out] = cqrrpt_mex(A, d_factor, eps, state) +// +// Outputs: +// Q m-by-k orthonormal factor, k = the rank CQRRPT detected +// R k-by-n upper-trapezoidal factor +// J 1-by-n vector of 1-based column pivots (int64), so +// A(:, J) ~= Q * R +// state_out RNG state after advancing +// +// k is discovered, not requested: CQRRPT is rank-revealing, and the rank it +// detects (its public `rank` member; eps acts only in its fallback +// re-estimation, see cqrrpt.m) sets the shared inner dimension of the +// returned Q and R. On a full-rank input k = n. +// +// J is 1-based, matching both MATLAB and CQRRPT's own convention: the +// pivots originate in LAPACK's geqp3 and RandLAPACK's util::col_swap reads +// them as `idx[i] - 1`. No conversion here. +// +// Shape contract worth noting while reading the buffer sizes below: CQRRPT +// factors A IN PLACE into Q, requiring an m-by-n work buffer even though only +// the first k columns come back as output, and it writes R into a caller- +// allocated n-by-n buffer (ldr >= n) of which only the first k rows are +// meaningful. Both outputs are therefore sliced out of oversized scratch, +// which is why neither can be factored directly into its output buffer the +// way bqrrp_mex's implicit mode does. + +#include "mex.hpp" +#include "mexAdapter.hpp" + +#include + +#include +#include +#include +#include +#include +#include + +using matlab::mex::ArgumentList; +using matlab::data::Array; +using matlab::data::ArrayFactory; +using matlab::data::ArrayType; +using matlab::data::StructArray; +template using TypedArray = matlab::data::TypedArray; + + +class MexFunction : public matlab::mex::Function { +private: + using RNG = RandBLAS::DefaultRNG; + static constexpr size_t CTR_LEN = 4; + static constexpr size_t KEY_LEN = 2; + + std::shared_ptr matlabPtr; + ArrayFactory factory; + + [[noreturn]] void raise(const std::string& id, const std::string& msg) { + matlabPtr->feval(u"error", 0, std::vector{ + factory.createCharArray(id), + factory.createCharArray(msg) + }); + std::terminate(); // feval(error) does not return; this is a backstop. + } + + RandBLAS::RNGState read_state(const Array& state_arr) { + if (state_arr.getType() != ArrayType::STRUCT) { + raise("randlapack:cqrrpt_mex:state", + "state must be a struct with .counter (uint32[4]) and .key (uint32[2])"); + } + StructArray s = state_arr; + if (s.getNumberOfElements() == 0) { + raise("randlapack:cqrrpt_mex:state", "state struct must be nonempty"); + } + Array ctr_field = s[0]["counter"]; + Array key_field = s[0]["key"]; + if (ctr_field.getType() != ArrayType::UINT32 + || ctr_field.getNumberOfElements() != CTR_LEN) { + raise("randlapack:cqrrpt_mex:state", + "state.counter must be a uint32 array of length 4"); + } + if (key_field.getType() != ArrayType::UINT32 + || key_field.getNumberOfElements() != KEY_LEN) { + raise("randlapack:cqrrpt_mex:state", + "state.key must be a uint32 array of length 2"); + } + TypedArray ctr = std::move(ctr_field); + TypedArray key = std::move(key_field); + RandBLAS::RNGState state; + for (size_t i = 0; i < CTR_LEN; ++i) state.counter.v[i] = ctr[i]; + for (size_t i = 0; i < KEY_LEN; ++i) state.key.v[i] = key[i]; + return state; + } + + Array write_state(const RandBLAS::RNGState& state) { + StructArray out = factory.createStructArray({1, 1}, {"counter", "key"}); + TypedArray ctr = factory.createArray({1, CTR_LEN}); + TypedArray key = factory.createArray({1, KEY_LEN}); + for (size_t i = 0; i < CTR_LEN; ++i) ctr[i] = state.counter.v[i]; + for (size_t i = 0; i < KEY_LEN; ++i) key[i] = state.key.v[i]; + out[0]["counter"] = std::move(ctr); + out[0]["key"] = std::move(key); + return out; + } + + template + void run_cqrrpt(ArgumentList& outputs, ArgumentList& inputs, + double d_factor_in, double eps_in) + { + const Array& A_in = inputs[0]; + auto dims = A_in.getDimensions(); + const int64_t m = static_cast(dims[0]); + const int64_t n = static_cast(dims[1]); + const size_t mn = static_cast(m) * static_cast(n); + const T d_factor = static_cast(d_factor_in); + const T eps = static_cast(eps_in); + + // Const view + one bulk memcpy: a non-const TypedArray would unshare + // (deep-copy) the input before any element is read. + const TypedArray A_typed = A_in; + const T* A_src = &*A_typed.cbegin(); + std::unique_ptr A_work(new T[mn]); + std::memcpy(A_work.get(), A_src, sizeof(T) * mn); + + // R is n-by-n scratch (ldr = n) regardless of the detected rank; J is + // written directly into its final output buffer. + std::vector R_work(static_cast(n) * static_cast(n)); + auto J_buf = factory.createBuffer(static_cast(n)); + + // J must be zero on entry: CQRRPT hands it to LAPACK's geqp3, which + // READS jpvt (a nonzero entry means "this column is fixed, move it to + // the front"), and createBuffer returns uninitialized memory. Current + // RandLAPACK zeroes J itself; zeroing here as well keeps the binding + // correct against installs that predate that, where the symptom was + // nondeterministic pivoting and a collapsing rank estimate. + std::fill(J_buf.get(), J_buf.get() + n, static_cast(0)); + + auto state = read_state(inputs[3]); + RandLAPACK::CQRRPT alg(false, eps); + // CQRRPT returns early on an all-zero input without setting `rank`, + // which its constructor leaves uninitialized. Start from 0 so that + // case reaches the zero-rank error below instead of a garbage size. + alg.rank = 0; + + int ret = alg.call(m, n, A_work.get(), m, R_work.data(), n, + J_buf.get(), d_factor, state); + if (ret != 0) { + raise("randlapack:cqrrpt_mex:returnCode", + "CQRRPT returned non-zero status " + std::to_string(ret)); + } + + const int64_t k = alg.rank; + if (k <= 0) { + raise("randlapack:cqrrpt_mex:zeroRank", + "CQRRPT detected rank 0: the sketch of A is exactly zero " + "(for example, A is all zeros)"); + } + + // Q = first k columns of the in-place factored A (column-major, so + // they are the leading m*k contiguous elements). + auto Q_buf = factory.createBuffer( + static_cast(m) * static_cast(k)); + std::memcpy(Q_buf.get(), A_work.get(), + sizeof(T) * static_cast(m) * static_cast(k)); + outputs[0] = factory.createArrayFromBuffer( + {static_cast(m), static_cast(k)}, std::move(Q_buf)); + + // R = first k rows of the n-by-n scratch, restrided from ldr=n to k. + auto R_buf = factory.createBuffer( + static_cast(k) * static_cast(n)); + T* R_out = R_buf.get(); + for (int64_t j = 0; j < n; ++j) { + std::memcpy(&R_out[j * k], &R_work[j * n], sizeof(T) * k); + } + outputs[1] = factory.createArrayFromBuffer( + {static_cast(k), static_cast(n)}, std::move(R_buf)); + + // J: 1-based int64 row vector, written in place by CQRRPT above. + outputs[2] = factory.createArrayFromBuffer( + {1, static_cast(n)}, std::move(J_buf)); + + if (outputs.size() >= 4) { + outputs[3] = write_state(state); + } + } + + double read_double(const Array& a, const char* name) { + if (a.getNumberOfElements() != 1) { + raise(std::string("randlapack:cqrrpt_mex:") + name, + std::string(name) + " must be a scalar"); + } + switch (a.getType()) { + case ArrayType::DOUBLE: return TypedArray(a)[0]; + case ArrayType::SINGLE: return static_cast(TypedArray(a)[0]); + case ArrayType::INT64: return static_cast(TypedArray(a)[0]); + default: + raise(std::string("randlapack:cqrrpt_mex:") + name, + std::string(name) + " must be a numeric scalar"); + } + } + +public: + MexFunction() : matlabPtr(getEngine()) {} + + void operator()(ArgumentList outputs, ArgumentList inputs) { + try { + if (inputs.size() != 4) { + raise("randlapack:cqrrpt_mex:nargin", + "Expected 4 input arguments: A, d_factor, eps, state"); + } + if (outputs.size() < 3 || outputs.size() > 4) { + raise("randlapack:cqrrpt_mex:nargout", + "Expected 3 or 4 output arguments"); + } + + const Array& A_in = inputs[0]; + if (A_in.getDimensions().size() != 2) { + raise("randlapack:cqrrpt_mex:A", + "A must be a 2D real numeric matrix (single or double)"); + } + const ArrayType A_type = A_in.getType(); + if (A_type != ArrayType::DOUBLE && A_type != ArrayType::SINGLE) { + raise("randlapack:cqrrpt_mex:type", "A must be single or double"); + } + + const double d_factor = read_double(inputs[1], "d_factor"); + const double eps = read_double(inputs[2], "eps"); + + if (d_factor < 1.0) { + raise("randlapack:cqrrpt_mex:d_factor", + "d_factor must be >= 1.0"); + } + if (eps < 0.0) { + raise("randlapack:cqrrpt_mex:eps", "eps must be nonnegative"); + } + + // CQRRPT sketches to d = d_factor*n rows and needs m >= d >= n. + // Checked here, with the arithmetic spelled out, rather than + // letting the requirement surface as a dimension complaint from + // three layers down inside the sketch-and-precondition chain. + const int64_t m = static_cast(A_in.getDimensions()[0]); + const int64_t n = static_cast(A_in.getDimensions()[1]); + const int64_t d = static_cast(d_factor * static_cast(n)); + if (m < d) { + raise("randlapack:cqrrpt_mex:notTall", + "CQRRPT requires a tall matrix: m >= d_factor*n. Got m=" + + std::to_string(m) + ", n=" + std::to_string(n) + + ", d_factor=" + std::to_string(d_factor) + + " (needs m >= " + std::to_string(d) + "). Reduce d_factor " + "toward 1.0, or use randlapack.bqrrp for wide matrices."); + } + + if (A_type == ArrayType::DOUBLE) { + run_cqrrpt(outputs, inputs, d_factor, eps); + } else { + run_cqrrpt(outputs, inputs, d_factor, eps); + } + } catch (const RandLAPACK::Error& e) { + raise("randlapack:cqrrpt_mex:RandLAPACKError", e.what()); + } catch (const RandBLAS::Error& e) { + raise("randlapack:cqrrpt_mex:RandBLASError", e.what()); + } catch (const matlab::engine::Exception&) { + // Errors raised by this MEX (raise() calls feval("error"), which + // throws matlab::engine::MATLABException) and other MATLAB engine + // errors propagate unchanged, so MATLAB sees the original error ID + // and message. In the MEX API (cppmex/mexException.hpp) these + // derive from std::exception but not from matlab::Exception (the + // C++ Engine API's classes of the same name differ): without this + // handler, the std::exception one below would re-raise them as + // ...:StdError. + throw; + } catch (const matlab::Exception&) { + // MATLAB Data API exceptions propagate unchanged. + throw; + } catch (const std::exception& e) { + raise("randlapack:cqrrpt_mex:StdError", e.what()); + } catch (...) { + raise("randlapack:cqrrpt_mex:Unknown", "Unknown exception in cqrrpt_mex"); + } + } +}; diff --git a/matlab/src/cqrrt_mex.cc b/matlab/src/cqrrt_mex.cc new file mode 100644 index 0000000..6bb789a --- /dev/null +++ b/matlab/src/cqrrt_mex.cc @@ -0,0 +1,260 @@ +// MEX wrapper for RandLAPACK::CQRRT using MATLAB's C++ Data API (R2018a+). +// +// MATLAB-side signature (low-level; users normally call randlapack.cqrrt.m): +// +// [Q, R, state_out] = cqrrt_mex(A, d_factor, compute_Q, state) +// +// Outputs: +// Q m-by-n orthonormal factor; 0-by-0 when compute_Q is false +// R n-by-n upper-triangular factor, so A ~= Q * R +// state_out RNG state after advancing +// +// CQRRT is unpivoted: it sketches A, takes R_sk from a QR of the sketch, +// preconditions A := A * inv(R_sk), and finishes with Cholesky QR. It expects +// a full-rank tall input and does not detect rank: a rank-deficient input +// usually completes, with R numerically singular. The driver returns a +// nonzero status only when the sketch's R has an exact zero on its diagonal +// (e.g. A has a zero column) or the Cholesky factorization fails; that is +// reported below as an error pointing at the rank-revealing alternatives. +// +// compute_Q = false selects CQRRT's R-only mode, which skips forming Q. The +// driver still overwrites A while preconditioning, so a work copy of the +// input is needed either way; with compute_Q = true that copy is the Q output +// buffer itself, and A is factored in place into it. + +#include "mex.hpp" +#include "mexAdapter.hpp" + +#include + +#include +#include +#include +#include +#include +#include + +using matlab::mex::ArgumentList; +using matlab::data::Array; +using matlab::data::ArrayFactory; +using matlab::data::ArrayType; +using matlab::data::StructArray; +template using TypedArray = matlab::data::TypedArray; + + +class MexFunction : public matlab::mex::Function { +private: + using RNG = RandBLAS::DefaultRNG; + static constexpr size_t CTR_LEN = 4; + static constexpr size_t KEY_LEN = 2; + + std::shared_ptr matlabPtr; + ArrayFactory factory; + + [[noreturn]] void raise(const std::string& id, const std::string& msg) { + matlabPtr->feval(u"error", 0, std::vector{ + factory.createCharArray(id), + factory.createCharArray(msg) + }); + std::terminate(); // feval(error) does not return; this is a backstop. + } + + RandBLAS::RNGState read_state(const Array& state_arr) { + if (state_arr.getType() != ArrayType::STRUCT) { + raise("randlapack:cqrrt_mex:state", + "state must be a struct with .counter (uint32[4]) and .key (uint32[2])"); + } + StructArray s = state_arr; + if (s.getNumberOfElements() == 0) { + raise("randlapack:cqrrt_mex:state", "state struct must be nonempty"); + } + Array ctr_field = s[0]["counter"]; + Array key_field = s[0]["key"]; + if (ctr_field.getType() != ArrayType::UINT32 + || ctr_field.getNumberOfElements() != CTR_LEN) { + raise("randlapack:cqrrt_mex:state", + "state.counter must be a uint32 array of length 4"); + } + if (key_field.getType() != ArrayType::UINT32 + || key_field.getNumberOfElements() != KEY_LEN) { + raise("randlapack:cqrrt_mex:state", + "state.key must be a uint32 array of length 2"); + } + TypedArray ctr = std::move(ctr_field); + TypedArray key = std::move(key_field); + RandBLAS::RNGState state; + for (size_t i = 0; i < CTR_LEN; ++i) state.counter.v[i] = ctr[i]; + for (size_t i = 0; i < KEY_LEN; ++i) state.key.v[i] = key[i]; + return state; + } + + Array write_state(const RandBLAS::RNGState& state) { + StructArray out = factory.createStructArray({1, 1}, {"counter", "key"}); + TypedArray ctr = factory.createArray({1, CTR_LEN}); + TypedArray key = factory.createArray({1, KEY_LEN}); + for (size_t i = 0; i < CTR_LEN; ++i) ctr[i] = state.counter.v[i]; + for (size_t i = 0; i < KEY_LEN; ++i) key[i] = state.key.v[i]; + out[0]["counter"] = std::move(ctr); + out[0]["key"] = std::move(key); + return out; + } + + template + void run_cqrrt(ArgumentList& outputs, ArgumentList& inputs, + double d_factor_in, bool compute_Q) + { + const Array& A_in = inputs[0]; + auto dims = A_in.getDimensions(); + const int64_t m = static_cast(dims[0]); + const int64_t n = static_cast(dims[1]); + const size_t mn = static_cast(m) * static_cast(n); + const size_t nn = static_cast(n) * static_cast(n); + const T d_factor = static_cast(d_factor_in); + + // Work copy of A: the Q output buffer when Q is wanted (factored in + // place, no second copy), otherwise scratch. Const view + one bulk + // memcpy: a non-const TypedArray would unshare (deep-copy) the input + // before any element is read. + const TypedArray A_typed = A_in; + const T* A_src = &*A_typed.cbegin(); + auto Q_buf = factory.createBuffer(compute_Q ? mn : 0); + std::unique_ptr A_scratch(compute_Q ? nullptr : new T[mn]); + T* A_work = compute_Q ? Q_buf.get() : A_scratch.get(); + std::memcpy(A_work, A_src, sizeof(T) * mn); + + // R must be zero on entry. CQRRT writes only the upper triangle until + // its final trmm, which multiplies the whole n-by-n buffer: anything + // left below the diagonal would leak into the upper triangle of R. + // createBuffer returns uninitialized memory, hence the fill. + auto R_buf = factory.createBuffer(nn); + std::fill(R_buf.get(), R_buf.get() + nn, static_cast(0)); + + auto state = read_state(inputs[3]); + // CQRRT keeps an eps member for interface symmetry with CQRRPT but + // does not read it. + RandLAPACK::CQRRT alg(false, static_cast(0)); + alg.compute_Q = compute_Q; + + int ret = alg.call(m, n, A_work, m, R_buf.get(), n, d_factor, state); + if (ret != 0) { + raise("randlapack:cqrrt_mex:breakdown", + "CQRRT could not factor A: the R factor of its sketch is " + "exactly singular (e.g. A has a zero column) or its Cholesky " + "step failed, so A is rank-deficient or too ill-conditioned " + "for an unpivoted factorization. Use randlapack.cqrrpt, which " + "reveals the rank, or randlapack.bqrrp."); + } + + if (compute_Q) { + outputs[0] = factory.createArrayFromBuffer( + {static_cast(m), static_cast(n)}, std::move(Q_buf)); + } else { + outputs[0] = factory.createArray({0, 0}); + } + if (outputs.size() >= 2) { + outputs[1] = factory.createArrayFromBuffer( + {static_cast(n), static_cast(n)}, std::move(R_buf)); + } + if (outputs.size() >= 3) { + outputs[2] = write_state(state); + } + } + + double read_double(const Array& a, const char* name) { + if (a.getNumberOfElements() != 1) { + raise(std::string("randlapack:cqrrt_mex:") + name, + std::string(name) + " must be a scalar"); + } + switch (a.getType()) { + case ArrayType::DOUBLE: return TypedArray(a)[0]; + case ArrayType::SINGLE: return static_cast(TypedArray(a)[0]); + case ArrayType::INT64: return static_cast(TypedArray(a)[0]); + default: + raise(std::string("randlapack:cqrrt_mex:") + name, + std::string(name) + " must be a numeric scalar"); + } + } + + bool read_logical(const Array& a, const char* name) { + if (a.getNumberOfElements() != 1 || a.getType() != ArrayType::LOGICAL) { + raise(std::string("randlapack:cqrrt_mex:") + name, + std::string(name) + " must be a logical scalar"); + } + return TypedArray(a)[0]; + } + +public: + MexFunction() : matlabPtr(getEngine()) {} + + void operator()(ArgumentList outputs, ArgumentList inputs) { + try { + if (inputs.size() != 4) { + raise("randlapack:cqrrt_mex:nargin", + "Expected 4 input arguments: A, d_factor, compute_Q, state"); + } + if (outputs.size() < 2 || outputs.size() > 3) { + raise("randlapack:cqrrt_mex:nargout", + "Expected 2 or 3 output arguments"); + } + + const Array& A_in = inputs[0]; + if (A_in.getDimensions().size() != 2) { + raise("randlapack:cqrrt_mex:A", + "A must be a 2D real numeric matrix (single or double)"); + } + const ArrayType A_type = A_in.getType(); + if (A_type != ArrayType::DOUBLE && A_type != ArrayType::SINGLE) { + raise("randlapack:cqrrt_mex:type", "A must be single or double"); + } + + const double d_factor = read_double(inputs[1], "d_factor"); + const bool compute_Q = read_logical(inputs[2], "compute_Q"); + if (d_factor < 1.0) { + raise("randlapack:cqrrt_mex:d_factor", "d_factor must be >= 1.0"); + } + + // CQRRT sketches to d = d_factor*n rows and needs m >= d >= n. + // Checked here, with the arithmetic spelled out, rather than + // letting the requirement surface as a dimension complaint from + // inside the sketch-and-precondition chain. + const int64_t m = static_cast(A_in.getDimensions()[0]); + const int64_t n = static_cast(A_in.getDimensions()[1]); + const int64_t d = static_cast(d_factor * static_cast(n)); + if (m < d) { + raise("randlapack:cqrrt_mex:notTall", + "CQRRT requires a tall matrix: m >= d_factor*n. Got m=" + + std::to_string(m) + ", n=" + std::to_string(n) + + ", d_factor=" + std::to_string(d_factor) + + " (needs m >= " + std::to_string(d) + "). Reduce d_factor " + "toward 1.0, or use randlapack.bqrrp for wide matrices."); + } + + if (A_type == ArrayType::DOUBLE) { + run_cqrrt(outputs, inputs, d_factor, compute_Q); + } else { + run_cqrrt(outputs, inputs, d_factor, compute_Q); + } + } catch (const RandLAPACK::Error& e) { + raise("randlapack:cqrrt_mex:RandLAPACKError", e.what()); + } catch (const RandBLAS::Error& e) { + raise("randlapack:cqrrt_mex:RandBLASError", e.what()); + } catch (const matlab::engine::Exception&) { + // Errors raised by this MEX (raise() calls feval("error"), which + // throws matlab::engine::MATLABException) and other MATLAB engine + // errors propagate unchanged, so MATLAB sees the original error ID + // and message. In the MEX API (cppmex/mexException.hpp) these + // derive from std::exception but not from matlab::Exception (the + // C++ Engine API's classes of the same name differ): without this + // handler, the std::exception one below would re-raise them as + // ...:StdError. + throw; + } catch (const matlab::Exception&) { + // MATLAB Data API exceptions propagate unchanged. + throw; + } catch (const std::exception& e) { + raise("randlapack:cqrrt_mex:StdError", e.what()); + } catch (...) { + raise("randlapack:cqrrt_mex:Unknown", "Unknown exception in cqrrt_mex"); + } + } +}; diff --git a/matlab/tests/test_cqrrpt.m b/matlab/tests/test_cqrrpt.m new file mode 100644 index 0000000..ea9fca5 --- /dev/null +++ b/matlab/tests/test_cqrrpt.m @@ -0,0 +1,181 @@ +function test_cqrrpt() +%test_cqrrpt Correctness checks for randlapack.cqrrpt. +% +% Runs CQRRPT on tall matrices and checks: +% 1. Orthogonality of Q. +% 2. Factorization fidelity: A(:, J) ~= Q*R. +% 3. R is upper-trapezoidal. +% 4. J is a valid 1-based permutation. +% 5. Rank detection on a deliberately rank-deficient input. +% 6. Reproducibility from an explicit RNG state. +% 7. Error paths, including the tall-matrix requirement. +% +% Run from MATLAB after: +% addpath('/path/to/randlapack-bindings/matlab') +% addpath('/path/to/randlapack-bindings/matlab/tests') +% test_cqrrpt + + types = {'double', 'single'}; + sizes = {[500, 50], [1000, 100], [2000, 200]}; + + for ti = 1:numel(types) + T = types{ti}; + for si = 1:numel(sizes) + sz = sizes{si}; + A = cast(randn(sz(1), sz(2)), T); + + [Q, R, J] = randlapack.cqrrpt(A); + + % Full-rank input: the detected rank should be n. + assert(size(R, 1) == sz(2), ... + '%s %dx%d: full-rank input gave rank %d, expected %d', ... + T, sz(1), sz(2), size(R, 1), sz(2)); + + check_orthogonality(Q, T); + check_factorization(A, Q, R, J, T); + check_upper_trapezoidal(R, T); + check_pivot_validity(J, sz(2)); + + fprintf('PASS: %s %dx%d full rank\n', T, sz(1), sz(2)); + end + + % --- Rank-deficient input: the rank must be revealed exactly --- + % Duplicated columns make the true rank exact and known, so this + % assertion is tight on purpose. It is also a regression test for a + % specific defect: CQRRPT passes J straight to LAPACK's geqp3, which + % READS jpvt on entry, so a J buffer that is not zeroed makes the + % pivoting depend on uninitialized memory. The visible symptom was + % exactly this check: the detected rank wandered between 1 and 40 + % across runs while the full-rank cases stayed green. + m = 2000; r = 40; extra = 10; + B = cast(randn(m, r), T); + A = [B, B(:, 1:extra)]; % rank r, r+extra columns + A = A(:, randperm(size(A, 2))); % scatter the duplicates + + [Q, R, J] = randlapack.cqrrpt(A); + detected = size(R, 1); + assert(detected == r, ... + '%s: expected rank %d on a rank-deficient input, got %d', ... + T, r, detected); + check_orthogonality(Q, T); + check_leading_factorization(A, Q, R, J, T); + % With the rank fully revealed, the whole permuted matrix + % reconstructs, not just the leading block. + check_factorization(A, Q, R, J, T); + fprintf('PASS: %s rank detection (%d of %d columns)\n', ... + T, detected, size(A, 2)); + end + + % --- Reproducibility from an explicit state --- + A = randn(1000, 100); + st = struct('counter', uint32([0 0 0 0]), 'key', uint32([0 11])); + [Q1, R1, J1] = randlapack.cqrrpt(A, 'state', st); + [Q2, R2, J2] = randlapack.cqrrpt(A, 'state', st); + assert(isequal(Q1, Q2) && isequal(R1, R2) && isequal(J1, J2), ... + 'same explicit state must give bit-identical results'); + fprintf('PASS: reproducible from an explicit RNG state\n'); + + [~, ~, ~, st_out] = randlapack.cqrrpt(A, 'state', st); + assert(~isequal(st_out.counter, st.counter) || ~isequal(st_out.key, st.key), ... + 'returned state must differ from the input state'); + fprintf('PASS: RNG state advances\n'); + + % --- eps is accepted (it only acts in the fallback rank estimate, so a + % full-rank input keeps its full rank at any eps) --- + [Q_e, R_e, J_e] = randlapack.cqrrpt(A, 'eps', 1e-6); + assert(size(R_e, 1) == 100, 'full-rank input must keep rank 100 with eps = 1e-6'); + check_factorization(A, Q_e, R_e, J_e, 'double'); + fprintf('PASS: eps option accepted\n'); + + % --- d_factor is honoured --- + [~, R_a, ~] = randlapack.cqrrpt(A, 'd_factor', 1.0); + [~, R_b, ~] = randlapack.cqrrpt(A, 'd_factor', 2.0); + assert(size(R_a, 1) == 100 && size(R_b, 1) == 100, ... + 'full-rank input must give rank 100 at either d_factor'); + fprintf('PASS: d_factor accepted across its range\n'); + + % --- Error paths --- + check_error(@() randlapack.cqrrpt(zeros(0, 5)), 'randlapack:cqrrpt:emptyInput'); + % An all-zero input has rank 0, which has no factorization to return. + check_error(@() randlapack.cqrrpt(zeros(500, 50)), 'randlapack:cqrrpt_mex:zeroRank'); + % Wide matrix: the tall requirement must be reported in MATLAB terms, + % not as a dimension complaint from inside the sketching chain. + check_error(@() randlapack.cqrrpt(randn(50, 500)), 'randlapack:cqrrpt:notTall'); + % Square is also too short once d_factor > 1. + check_error(@() randlapack.cqrrpt(randn(100, 100)), 'randlapack:cqrrpt:notTall'); + check_error(@() randlapack.cqrrpt(A, 'd_factor', 0.5), ... + 'MATLAB:cqrrpt:notGreaterEqual'); + check_error(@() randlapack.cqrrpt(A, 'nosuchopt', 1), ... + 'randlapack:cqrrpt:unknownOption'); + check_error(@() randlapack.cqrrpt(int32(A)), 'MATLAB:cqrrpt:invalidType'); + fprintf('PASS: error paths\n'); + + fprintf('test_cqrrpt: all checks passed\n'); +end + + +function check_orthogonality(Q, T) + [m, k] = size(Q); + I = cast(eye(k), T); + err = norm(double(Q)' * double(Q) - double(I), 'fro'); + tol = sqrt(double(m * k)) * eps(T) * 100; + assert(err < tol, ... + 'orthogonality check failed: ||Q^T Q - I||_F = %g > tol = %g', ... + err, tol); +end + + +function check_factorization(A, Q, R, J, T) + [m, n] = size(A); + err = norm(double(A(:, J)) - double(Q) * double(R), 'fro') ... + / norm(double(A), 'fro'); + tol = sqrt(double(m * n)) * eps(T) * 100; + assert(err < tol, ... + 'factorization check failed: ||A(:,J) - QR||_F/||A||_F = %g > tol = %g', ... + err, tol); +end + + +function check_leading_factorization(A, Q, R, J, T) +% The leading k pivoted columns are reproduced exactly whatever rank came +% back, because R(:, 1:k) is precisely their triangular factor. This is the +% invariant that survives a conservative rank estimate. + [m, ~] = size(A); + k = size(R, 1); + lead = double(A(:, J(1:k))); + err = norm(lead - double(Q) * double(R(:, 1:k)), 'fro') ... + / norm(double(A), 'fro'); + tol = sqrt(double(m * k)) * eps(T) * 100; + assert(err < tol, ... + ['leading-block factorization failed: ' ... + '||A(:,J(1:k)) - Q*R(:,1:k)||_F/||A||_F = %g > tol = %g'], ... + err, tol); +end + + +function check_upper_trapezoidal(R, T) + below = tril(double(R), -1); + assert(norm(below, 'fro') == 0, ... + 'R must be upper-trapezoidal; found %g below the diagonal', ... + norm(below, 'fro')); + assert(isa(R, T), 'R must have the same class as A'); +end + + +function check_pivot_validity(J, n) + assert(numel(J) == n, 'J must have n entries'); + assert(isequal(sort(double(J(:)))', 1:n), ... + 'J must be a permutation of 1:n (1-based)'); +end + + +function check_error(fn, expected_id) + try + fn(); + catch err + assert(strcmp(err.identifier, expected_id), ... + 'expected error %s, got %s', expected_id, err.identifier); + return; + end + error('expected error %s, but no error was raised', expected_id); +end diff --git a/matlab/tests/test_cqrrt.m b/matlab/tests/test_cqrrt.m new file mode 100644 index 0000000..f737fed --- /dev/null +++ b/matlab/tests/test_cqrrt.m @@ -0,0 +1,177 @@ +function test_cqrrt() +%test_cqrrt Correctness checks for randlapack.cqrrt. +% +% Runs CQRRT on tall full-rank matrices and checks: +% 1. Orthogonality of Q. +% 2. Factorization fidelity: A ~= Q*R. +% 3. R is upper-triangular. +% 4. The R-only call returns the same R as the [Q, R] call. +% 5. An ill-conditioned (but full-rank) input still factors accurately. +% 6. An exactly singular sketch (zero column) is reported as an error; +% a merely rank-deficient input completes with a singular R. +% 7. Reproducibility from an explicit RNG state. +% 8. Error paths, including the tall-matrix requirement. +% +% Run from MATLAB after: +% addpath('/path/to/randlapack-bindings/matlab') +% addpath('/path/to/randlapack-bindings/matlab/tests') +% test_cqrrt + + types = {'double', 'single'}; + sizes = {[500, 50], [1000, 100], [2000, 200]}; + + for ti = 1:numel(types) + T = types{ti}; + for si = 1:numel(sizes) + sz = sizes{si}; + A = cast(randn(sz(1), sz(2)), T); + + [Q, R] = randlapack.cqrrt(A); + + assert(isequal(size(Q), sz) && isequal(size(R), [sz(2), sz(2)]), ... + '%s %dx%d: Q is %dx%d and R is %dx%d', T, sz(1), sz(2), ... + size(Q, 1), size(Q, 2), size(R, 1), size(R, 2)); + check_orthogonality(Q, T); + check_factorization(A, Q, R, T); + check_upper_triangular(R, T); + + fprintf('PASS: %s %dx%d\n', T, sz(1), sz(2)); + end + + % --- R-only mode returns the same R --- + % R does not depend on whether Q is formed: compute_Q only switches + % the triangular solve that forms Q, which reads the Cholesky factor + % but never writes it. From the same RNG state the two calls agree + % up to rounding. Not bit for bit: the two calls run + % on different work buffers, and an optimized BLAS may pick a + % different code path for a different memory alignment. + A = cast(randn(1000, 100), T); + st = struct('counter', uint32([0 0 0 0]), 'key', uint32([0 7])); + R_only = randlapack.cqrrt(A, 'state', st); + [~, R_full] = randlapack.cqrrt(A, 'state', st); + diff_R = norm(double(R_only) - double(R_full), 'fro') ... + / norm(double(R_full), 'fro'); + assert(diff_R <= 100 * eps(T), ... + '%s: R-only and [Q, R] calls differ by %g relative', T, diff_R); + assert(isa(R_only, T), 'R must have the same class as A'); + fprintf('PASS: %s R-only mode matches [Q, R]\n', T); + + % --- An exactly singular sketch is reported as an error --- + % A zero column makes the sketch's R exactly singular, so CQRRT + % cannot precondition and returns a failure status. + A = cast(randn(1000, 100), T); + A(:, 37) = 0; + check_error(@() randlapack.cqrrt(A), 'randlapack:cqrrt_mex:breakdown'); + fprintf('PASS: %s error on an exactly singular sketch (zero column)\n', T); + end + + % --- A rank-deficient input is NOT reported (documented behaviour) --- + % Duplicated columns: the sketch's R diagonal is tiny but not zero, so + % CQRRT completes. Q*R still reproduces A, but R is numerically + % singular, which is what the help tells users to check for. + B = randn(2000, 40); + A = [B, B(:, 1:10)]; + A = A(:, randperm(50)); + [Q, R] = randlapack.cqrrt(A); + d = abs(diag(R)); + assert(min(d) / max(d) < 1e-10, ... + 'rank-deficient input: expected a numerically singular R, got min/max diag %g', ... + min(d) / max(d)); + check_factorization(A, Q, R, 'double'); + fprintf('PASS: double rank-deficient input completes with a singular R\n'); + + % --- Ill-conditioned but full-rank input (double) --- + % Sketch-based preconditioning is what lets Cholesky QR survive here: + % plain Cholesky QR on A'A would square cond(A) = 1e10 past 1/eps. + m = 4000; n = 100; + [U, ~] = qr(randn(m, n), 0); + [V, ~] = qr(randn(n, n)); + A = U * diag(logspace(0, -10, n)) * V'; + [Q, R] = randlapack.cqrrt(A); + check_orthogonality(Q, 'double'); + check_factorization(A, Q, R, 'double'); + fprintf('PASS: double, cond(A) = 1e10\n'); + + % --- Reproducibility from an explicit state --- + A = randn(1000, 100); + st = struct('counter', uint32([0 0 0 0]), 'key', uint32([0 11])); + [Q1, R1] = randlapack.cqrrt(A, 'state', st); + [Q2, R2] = randlapack.cqrrt(A, 'state', st); + assert(isequal(Q1, Q2) && isequal(R1, R2), ... + 'same explicit state must give bit-identical results'); + fprintf('PASS: reproducible from an explicit RNG state\n'); + + [~, ~, st_out] = randlapack.cqrrt(A, 'state', st); + assert(~isequal(st_out.counter, st.counter) || ~isequal(st_out.key, st.key), ... + 'returned state must differ from the input state'); + fprintf('PASS: RNG state advances\n'); + + % --- d_factor is honoured --- + [Q_a, R_a] = randlapack.cqrrt(A, 'd_factor', 1.0); + [Q_b, R_b] = randlapack.cqrrt(A, 'd_factor', 2.0); + check_factorization(A, Q_a, R_a, 'double'); + check_factorization(A, Q_b, R_b, 'double'); + fprintf('PASS: d_factor accepted across its range\n'); + + % --- Error paths --- + check_error(@() randlapack.cqrrt(zeros(0, 5)), 'randlapack:cqrrt:emptyInput'); + % Wide matrix: the tall requirement must be reported in MATLAB terms, + % not as a dimension complaint from inside the sketching chain. + check_error(@() randlapack.cqrrt(randn(50, 500)), 'randlapack:cqrrt:notTall'); + % Square is also too short once d_factor > 1. + check_error(@() randlapack.cqrrt(randn(100, 100)), 'randlapack:cqrrt:notTall'); + check_error(@() randlapack.cqrrt(A, 'd_factor', 0.5), ... + 'MATLAB:cqrrt:notGreaterEqual'); + check_error(@() randlapack.cqrrt(A, 'nosuchopt', 1), ... + 'randlapack:cqrrt:unknownOption'); + % CQRRT has no eps option (the driver does not use one). + check_error(@() randlapack.cqrrt(A, 'eps', 1e-8), ... + 'randlapack:cqrrt:unknownOption'); + check_error(@() randlapack.cqrrt(int32(A)), 'MATLAB:cqrrt:invalidType'); + fprintf('PASS: error paths\n'); + + fprintf('test_cqrrt: all checks passed\n'); +end + + +function check_orthogonality(Q, T) + [m, n] = size(Q); + I = cast(eye(n), T); + err = norm(double(Q)' * double(Q) - double(I), 'fro'); + tol = sqrt(double(m * n)) * eps(T) * 100; + assert(err < tol, ... + 'orthogonality check failed: ||Q^T Q - I||_F = %g > tol = %g', ... + err, tol); +end + + +function check_factorization(A, Q, R, T) + [m, n] = size(A); + err = norm(double(A) - double(Q) * double(R), 'fro') ... + / norm(double(A), 'fro'); + tol = sqrt(double(m * n)) * eps(T) * 100; + assert(err < tol, ... + 'factorization check failed: ||A - QR||_F/||A||_F = %g > tol = %g', ... + err, tol); +end + + +function check_upper_triangular(R, T) + below = tril(double(R), -1); + assert(norm(below, 'fro') == 0, ... + 'R must be upper-triangular; found %g below the diagonal', ... + norm(below, 'fro')); + assert(isa(R, T), 'R must have the same class as A'); +end + + +function check_error(fn, expected_id) + try + fn(); + catch err + assert(strcmp(err.identifier, expected_id), ... + 'expected error %s, got %s', expected_id, err.identifier); + return; + end + error('expected error %s, but no error was raised', expected_id); +end