Require a square matrix in solve(), matching R

solve(a, b) with a rectangular a fell through to a least-squares dgels
call, returning qr.solve()'s answer where R raises "'a' (m x n) must be
square". Statically rectangular systems are now a compile error and
symbolic squareness is guarded at run time before the dgesv call, via
the assert_square_matrix() helper the other LAPACK lowerings already
use. The now-unreachable rectangular tail of lapack_solve() (dgels, and
a dgelsy branch qr.solve() never reached) is deleted; qr.solve() keeps
its least-squares behavior.

The solve output follows ncol(a) while b follows nrow(a); when ncol is
statically 1 the output declares as a Fortran scalar, so a
symbolic-length b is copied elementwise instead of by whole-array
assignment.

Author: mns-nordicals

NEWS.md

@@ -32,6 +32,13 @@
   build even in arithmetic, and a symbolic-length `x + m` returned a plain
   vector where R returns a `1x1` matrix.
 
+- `solve(a, b)` now requires a square `a`, matching R. A rectangular `a`
+  previously fell through to a least-squares solve (dgels), returning an
+  answer where R raises `'a' (m x n) must be square`. Statically
+  rectangular systems are now a compile-time error; when squareness is not
+  known at compile time it is checked at run time. Use `qr.solve()` for
+  least-squares solutions of rectangular systems (unchanged).
+
 # quickr 0.3.0
 
 This release adds major new support for linear algebra, local functions,

R/r2f-matrix-blas.R

@@ -643,13 +643,12 @@ lapack_solve <- function(
 
   nrhs <- if (b_rank == 1L) 1L else dim_or_one(B, 2L)
 
-  # solve(a, b) with a rectangular `a` deliberately falls through to the
-  # least-squares branch below -- a divergence from base R (which requires
-  # a square `a`), locked by the "least-squares" tests in
-  # test-matrix-lapack.R. Squareness is a routing decision here, not a
-  # correctness guard: unknown squareness routes to dgels, which solves
-  # square systems exactly too.
-  if (dims_match(m, n) && !identical(context, "qr.solve")) {
+  # R's solve() requires a square `a`; least squares is qr.solve()'s job.
+  # Statically rectangular `a` is a compile error, symbolic dims get a
+  # runtime guard before the dgesv call. (A rectangular `a` used to fall
+  # through to a dgels least-squares solve -- an answer where R errors.)
+  if (!identical(context, "qr.solve")) {
+    assert_square_matrix(a_dims, A, context, hoist, scope)
     A_work <- hoist$declare_tmp(mode = "double", dims = list(m, m))
     hoist$emit(glue("{A_work@name} = {A_name}"))
 
@@ -671,7 +670,17 @@ lapack_solve <- function(
       out_var <- hoist$declare_tmp(mode = "double", dims = expected_dims)
       out_name <- out_var@name
     }
-    hoist$emit(glue("{out_name} = {B_input_name}"))
+    # The output length follows ncol(a) (R's contract) while `b` follows
+    # nrow(a); the two are only runtime-equal. When ncol is statically 1
+    # the output declares as a scalar, so a symbolic-length `b` must be
+    # copied elementwise, not by whole-array assignment.
+    b_src <- if (passes_as_scalar(out_var) && !passes_as_scalar(B@value)) {
+      subs <- str_flatten_commas(rep("1", b_rank))
+      glue("{B_input_name}({subs})")
+    } else {
+      B_input_name
+    }
+    hoist$emit(glue("{out_name} = {b_src}"))
 
     ipiv <- hoist$declare_tmp(mode = "integer", dims = list(m))
     info <- hoist$declare_tmp(mode = "integer", dims = NULL)
@@ -813,121 +822,6 @@ end do"
     }
     return(out)
   }
-
-  A_work <- hoist$declare_tmp(mode = "double", dims = list(m, n))
-  hoist$emit(glue("{A_work@name} = {A_name}"))
-
-  max_mn <- call("max", m, n)
-
-  B_work <- hoist$declare_tmp(mode = "double", dims = list(max_mn, nrhs))
-  m_f <- dims2f(list(m), scope)
-  if (!nzchar(m_f)) {
-    m_f <- "1"
-  }
-  n_f <- dims2f(list(n), scope)
-  if (!nzchar(n_f)) {
-    n_f <- "1"
-  }
-  nrhs_f <- dims2f(list(nrhs), scope)
-  if (!nzchar(nrhs_f)) {
-    nrhs_f <- "1"
-  }
-  hoist$emit(glue("{B_work@name} = 0.0_c_double"))
-  if (b_rank == 1L) {
-    hoist$emit(glue("{B_work@name}(1:{m_f}, 1) = {B_input_name}"))
-  } else {
-    hoist$emit(glue("{B_work@name}(1:{m_f}, 1:{nrhs_f}) = {B_input_name}"))
-  }
-
-  info <- hoist$declare_tmp(mode = "integer", dims = NULL)
-
-  mn <- call("min", m, n)
-  if (identical(context, "qr.solve")) {
-    jpvt <- hoist$declare_tmp(mode = "integer", dims = list(n))
-    hoist$emit(glue("{jpvt@name} = 0_c_int"))
-
-    rcond <- if (is.null(tol)) "1e-7_c_double" else as.character(tol)
-    rank <- hoist$declare_tmp(mode = "integer", dims = NULL)
-
-    lwork <- call(
-      "max",
-      1L,
-      call("+", mn, call("max", mn, nrhs)),
-      call("+", call("*", 2L, mn), call("*", 64L, call("+", n, 1L))),
-      call("+", mn, call("*", 2L, n))
-    )
-    work <- hoist$declare_tmp(mode = "double", dims = list(lwork))
-
-    hoist$emit(glue(
-      "call dgelsy({blas_int(m)}, {blas_int(n)}, {blas_int(nrhs)}, {A_work@name}, {blas_int(m)}, {B_work@name}, {blas_int(max_mn)}, {jpvt@name}, {rcond}, {rank@name}, {work@name}, {blas_int(lwork)}, {info@name})"
-    ))
-    emit_quickr_error_if(
-      condition = glue("{info@name} < 0_c_int"),
-      message = "Lapack routine dgelsy: illegal argument",
-      hoist = hoist,
-      scope = scope
-    )
-    emit_quickr_error_if(
-      condition = glue("{info@name} > 0_c_int"),
-      message = "Lapack routine dgelsy failed to converge",
-      hoist = hoist,
-      scope = scope
-    )
-    emit_quickr_error_if(
-      condition = glue("{rank@name} < {blas_int(n)}"),
-      message = "rank deficient matrix in qr.solve",
-      hoist = hoist,
-      scope = scope
-    )
-  } else {
-    lwork <- call("max", 1L, call("+", mn, call("max", mn, nrhs)))
-    work <- hoist$declare_tmp(mode = "double", dims = list(lwork))
-
-    hoist$emit(glue(
-      "call dgels('N', {blas_int(m)}, {blas_int(n)}, {blas_int(nrhs)}, {A_work@name}, {blas_int(m)}, {B_work@name}, {blas_int(max_mn)}, {work@name}, {blas_int(lwork)}, {info@name})"
-    ))
-    emit_quickr_error_if(
-      condition = glue("{info@name} < 0_c_int"),
-      message = "Lapack routine dgels: illegal argument",
-      hoist = hoist,
-      scope = scope
-    )
-  }
-
-  expected_dims <- if (b_rank == 1L) list(n) else list(n, nrhs)
-  writes_to_dest <- FALSE
-  if (
-    can_use_output(
-      dest,
-      input_names = c(A_name, B_input_name),
-      expected_dims = expected_dims,
-      context = context,
-      allow_alias = B_input_name
-    )
-  ) {
-    out_var <- dest
-    out_name <- dest@name
-    writes_to_dest <- TRUE
-  } else {
-    out_var <- hoist$declare_tmp(mode = "double", dims = expected_dims)
-    out_name <- out_var@name
-  }
-
-  if (b_rank == 1L) {
-    if (passes_as_scalar(out_var)) {
-      hoist$emit(glue("{out_name} = {B_work@name}(1, 1)"))
-    } else {
-      hoist$emit(glue("{out_name} = {B_work@name}(1:{n_f}, 1)"))
-    }
-  } else {
-    hoist$emit(glue("{out_name} = {B_work@name}(1:{n_f}, 1:{nrhs_f})"))
-  }
-
-  out <- Fortran(out_name, out_var)
-  if (writes_to_dest) {
-    out@writes_to_dest <- TRUE
-  }
-  out
 }
 
 lapack_inverse <- function(A, scope, hoist, dest = NULL, context = "solve") {

tests/testthat/test-matrix-lapack.R

@@ -62,66 +62,55 @@ test_that("solve handles column RHS matrices and 1x1 systems", {
   expect_quick_equal(solve_scalar, list(A = matrix(2.5, 1L, 1L), b = 1.25))
 })
 
-test_that("solve supports least-squares for rectangular systems", {
-  solve_ls_vec <- function(X, y) {
+test_that("solve requires a square coefficient matrix, like R", {
+  # Squareness unknown at compile time: a runtime guard runs before dgesv.
+  solve_sym <- function(X, y) {
     declare(
       type(X = double(n, k)),
       type(y = double(n))
     )
     solve(X, y)
   }
 
-  solve_ls_mat <- function(X, Y) {
-    declare(
-      type(X = double(n, k)),
-      type(Y = double(n, p))
-    )
-    solve(X, Y)
-  }
-
   set.seed(123)
   n <- 20
   k <- 5
-  p <- 3
   X <- matrix(rnorm(n * k), n, k)
   y <- rnorm(n)
-  Y <- matrix(rnorm(n * p), n, p)
 
-  q_solve_ls_vec <- expect_warning(quick(solve_ls_vec), NA)
-  q_solve_ls_mat <- expect_warning(quick(solve_ls_mat), NA)
+  q_solve_sym <- quick(solve_sym)
+  expect_error(q_solve_sym(X, y), "solve requires a square matrix")
+  expect_error(solve(X, y), "must be square") # the R oracle errors too
 
-  expect_equal(q_solve_ls_vec(X, y), qr.solve(X, y))
-  expect_equal(q_solve_ls_mat(X, Y), qr.solve(X, Y))
-})
+  base <- matrix(rnorm(n * n), n, n)
+  A <- crossprod(base) + diag(n)
+  expect_equal(q_solve_sym(A, y), solve(A, y))
 
-test_that("solve supports least-squares for single-column systems", {
-  solve_ls_col <- function(X, y) {
+  # One statically known axis still guards the symbolic one.
+  solve_col <- function(X, y) {
     declare(
       type(X = double(n, 1L)),
       type(y = double(n))
     )
     solve(X, y)
   }
-
-  set.seed(125)
-  n <- 20
-  X <- matrix(rnorm(n), n, 1L)
-  y <- rnorm(n)
-
-  q_solve_ls_col <- expect_warning(quick(solve_ls_col), NA)
-  expect_equal(q_solve_ls_col(X, y), qr.solve(X, y))
+  q_solve_col <- quick(solve_col)
+  expect_error(q_solve_col(X[, 1L, drop = FALSE], y), "square matrix")
+  expect_equal(
+    q_solve_col(matrix(2.5, 1L, 1L), 1.25),
+    solve(matrix(2.5, 1L, 1L), 1.25)
+  )
 })
 
-test_that("solve compiles 1-row least-squares systems", {
+test_that("solve rejects statically rectangular systems at compile time", {
   solve_one_row <- function(X, y) {
     declare(
       type(X = double(1L, 2L)),
       type(y = double(1L))
     )
     solve(X, y)
   }
-
-  expect_no_error(r2f(solve_one_row))
+  expect_error(r2f(solve_one_row), "solve requires a square matrix")
 })
 
 test_that("qr.solve matches R for vectors and matrices", {
@@ -195,7 +184,7 @@ test_that("qr.solve uses QR with pivoting for known square systems", {
   expect_false(has_call(square_fortran, "dgesv"))
 })
 
-test_that("solve uses dgesv for known square and dgels for rectangular systems", {
+test_that("solve always uses dgesv, guarding squareness when symbolic", {
   solve_square <- function(A, b) {
     declare(
       type(A = double(n, n)),
@@ -212,12 +201,7 @@ test_that("solve uses dgesv for known square and dgels for rectangular systems",
     solve(A, b)
   }
 
-  solve_rect <- function(A, b) {
-    declare(type(A = double(3, 2)), type(b = double(3)))
-    solve(A, b)
-  }
-
-  solve_rect_named <- function(A, b) {
+  solve_sym <- function(A, b) {
     declare(
       type(A = double(n, k)),
       type(b = double(n))
@@ -233,32 +217,27 @@ test_that("solve uses dgesv for known square and dgels for rectangular systems",
     capture.output(cat(r2f(solve_square_fixed))),
     collapse = "\n"
   )
-  rect_fixed_fortran <- paste(
-    capture.output(cat(r2f(solve_rect))),
-    collapse = "\n"
-  )
-  rect_named_fortran <- paste(
-    capture.output(cat(r2f(solve_rect_named))),
+  sym_fortran <- paste(
+    capture.output(cat(r2f(solve_sym))),
     collapse = "\n"
   )
 
   has_call <- function(code, routine) {
     any(grepl(paste0("call ", routine, "("), tolower(code), fixed = TRUE))
   }
 
-  # Search emitted Fortran to ensure solve() chooses LU (dgesv) for proven-square
-  # systems and least-squares (dgels) for rectangular systems.
+  # solve() is LU (dgesv) only; least squares is qr.solve()'s job. Provably
+  # square systems get no guard, symbolic squareness is checked at run time.
   expect_true(has_call(square_named_fortran, "dgesv"))
   expect_false(has_call(square_named_fortran, "dgels"))
+  expect_false(grepl("square matrix", square_named_fortran, fixed = TRUE))
 
   expect_true(has_call(square_fixed_fortran, "dgesv"))
   expect_false(has_call(square_fixed_fortran, "dgels"))
 
-  expect_true(has_call(rect_fixed_fortran, "dgels"))
-  expect_false(has_call(rect_fixed_fortran, "dgesv"))
-
-  expect_true(has_call(rect_named_fortran, "dgels"))
-  expect_false(has_call(rect_named_fortran, "dgesv"))
+  expect_true(has_call(sym_fortran, "dgesv"))
+  expect_false(has_call(sym_fortran, "dgels"))
+  expect_match(sym_fortran, "solve requires a square matrix", fixed = TRUE)
 })
 
 test_that("chol and chol2inv match R", {
@@ -415,7 +394,7 @@ test_that("Test bad path in lapack functions", {
   }
 
   expect_error(quick(solve_bad_rank), "expects a matrix")
-  expect_warning(quick(solve_non_square), NA)
+  expect_error(quick(solve_non_square), "requires a square matrix")
   expect_error(quick(solve_bad_rhs), "only supports vector or matrix")
   expect_error(quick(chol_bad_rank), "expects a matrix")
   expect_error(quick(chol_pivot), "pivot = TRUE")
@@ -446,7 +425,7 @@ test_that("lapack functions test non-square matrix errors", {
     chol2inv(R)
   }
 
-  expect_warning(quick(solve_rect), NA)
+  expect_error(quick(solve_rect), "requires a square matrix")
   expect_error(quick(chol_rect), "requires a square matrix")
   expect_error(quick(chol2inv_rect), "requires a square matrix")
 })