Refuse complex operands in linear algebra

The BLAS/LAPACK paths cast logical/integer operands to double via
maybe_cast_double(), which passes complex through untouched -- so
complex operands flowed into the real d* routines, which read complex
storage as reals: complex(2) %*% complex(2) returned a real dot product
of the real parts where R returns the complex result. A silent wrong
answer, found in external review.

All linalg operand casts now go through cast_linalg_double(), which
refuses complex with a compile-time message naming the divergence
(linear algebra in quickr is double-only). Covers %*% (including t()
forms), crossprod/tcrossprod, solve/qr.solve, forwardsolve/backsolve,
chol, chol2inv, svd, and outer. The mode-preserving standalone t() and
elementwise complex arithmetic are untouched.

Author: mns-nordicals

R/r2f-matrix-blas.R

@@ -471,8 +471,8 @@ outer_mul <- function(
 ) {
   assert_hoist_env(hoist)
 
-  x <- maybe_cast_double(x)
-  y <- maybe_cast_double(y)
+  x <- cast_linalg_double(x, context)
+  y <- cast_linalg_double(y, context)
 
   if (x@value@rank > 1L || y@value@rank > 1L) {
     stop("outer() only supports vectors or scalars")
@@ -523,8 +523,8 @@ triangular_solve <- function(
 ) {
   assert_hoist_env(hoist)
 
-  A <- maybe_cast_double(A)
-  B <- maybe_cast_double(B)
+  A <- cast_linalg_double(A, context)
+  B <- cast_linalg_double(B, context)
 
   assert_rank2_matrix(A, "triangular solve expects a matrix")
 
@@ -605,8 +605,8 @@ lapack_solve <- function(
 ) {
   assert_hoist_env(hoist)
 
-  A <- maybe_cast_double(A)
-  B <- maybe_cast_double(B)
+  A <- cast_linalg_double(A, context)
+  B <- cast_linalg_double(B, context)
 
   assert_rank2_matrix(A, paste0(context, " expects a matrix for `a`"))
 
@@ -812,7 +812,7 @@ end do"
 lapack_inverse <- function(A, scope, hoist, dest = NULL, context = "solve") {
   assert_hoist_env(hoist)
 
-  A <- maybe_cast_double(A)
+  A <- cast_linalg_double(A, context)
   assert_rank2_matrix(A, paste0(context, " expects a matrix for `a`"))
 
   a_dims <- matrix_dims(A)
@@ -886,7 +886,7 @@ lapack_inverse <- function(A, scope, hoist, dest = NULL, context = "solve") {
 lapack_chol <- function(A, scope, hoist, dest = NULL, context = "chol") {
   assert_hoist_env(hoist)
 
-  A <- maybe_cast_double(A)
+  A <- cast_linalg_double(A, context)
   assert_rank2_matrix(A, paste0(context, " expects a matrix"))
 
   a_dims <- matrix_dims(A)
@@ -949,7 +949,7 @@ lapack_chol2inv <- function(
 ) {
   assert_hoist_env(hoist)
 
-  R <- maybe_cast_double(R)
+  R <- cast_linalg_double(R, context)
   assert_rank2_matrix(R, paste0(context, " expects a matrix"))
 
   r_dims <- matrix_dims(R)
@@ -1161,7 +1161,7 @@ lapack_svd <- function(
   assert_hoist_env(hoist)
   stopifnot(inherits(d, Variable), inherits(u, Variable), inherits(v, Variable))
 
-  A <- maybe_cast_double(A)
+  A <- cast_linalg_double(A, context)
   dims <- svd_dims(A, context = context)
   m <- dims$m
   n <- dims$n

R/r2f-matrix-parse.R

@@ -9,7 +9,7 @@ unwrap_transpose_arg <- function(arg, scope, ..., hoist) {
   if (is_call(arg_unwrapped, quote(t)) && length(arg_unwrapped) == 2L) {
     inner_arg <- unwrap_parens(arg_unwrapped[[2L]])
     inner <- r2f(inner_arg, scope, ..., hoist = hoist)
-    inner <- maybe_cast_double(inner)
+    inner <- cast_linalg_double(inner, "%*%")
     if (inner@value@rank == 2L) {
       return(list(value = inner, trans = "T"))
     } else if (inner@value@rank == 1L) {
@@ -26,7 +26,7 @@ unwrap_transpose_arg <- function(arg, scope, ..., hoist) {
     }
   }
   value <- r2f(arg, scope, ..., hoist = hoist)
-  value <- maybe_cast_double(value)
+  value <- cast_linalg_double(value, "%*%")
   list(value = value, trans = "N")
 }
 

R/r2f-matrix.R

@@ -627,7 +627,10 @@ register_r2f_handler(
     tol <- if (is.null(tol_arg) || is_missing(tol_arg)) {
       r2f(1e-7, scope, ..., hoist = hoist)
     } else {
-      tol <- maybe_cast_double(r2f(tol_arg, scope, ..., hoist = hoist))
+      tol <- cast_linalg_double(
+        r2f(tol_arg, scope, ..., hoist = hoist),
+        "qr.solve"
+      )
       if (tol@value@rank != 0L) {
         stop("qr.solve() expects a scalar `tol`", call. = FALSE)
       }
@@ -916,7 +919,7 @@ crossprod_like <- function(
   context
 ) {
   x <- r2f(x_arg, scope, ..., hoist = hoist)
-  x <- maybe_cast_double(x)
+  x <- cast_linalg_double(x, context)
 
   if (is.null(y_arg)) {
     return(syrk(
@@ -929,7 +932,7 @@ crossprod_like <- function(
     ))
   }
 
-  y <- maybe_cast_double(r2f(y_arg, scope, ..., hoist = hoist))
+  y <- cast_linalg_double(r2f(y_arg, scope, ..., hoist = hoist), context)
 
   x_dims <- matrix_dims(x)
   y_dims <- matrix_dims(y)

R/r2f-operators-helpers.R

@@ -85,6 +85,24 @@ maybe_cast_double <- function(x) {
   }
 }
 
+# Cast a linear-algebra operand to double for the real BLAS/LAPACK
+# lowerings (dgemm, dgesv, ...). Complex operands are refused: the d*
+# routines would read complex storage as reals and return a plausible
+# wrong answer, and quickr has no z* lowerings. R supports complex
+# linear algebra, so the message names the divergence.
+# Used by: r2f-matrix.R, r2f-matrix-parse.R, r2f-matrix-blas.R
+cast_linalg_double <- function(x, context) {
+  if (identical(x@value@mode, "complex")) {
+    stop(
+      context,
+      " does not support complex operands; ",
+      "linear algebra in quickr is double-only",
+      call. = FALSE
+    )
+  }
+  maybe_cast_double(x)
+}
+
 # Promote a list of operands to their common (lattice-join) mode, casting
 # each one whose mode differs. For contexts where Fortran requires uniform
 # argument types: array constructors (c()), min/max, merge, modulo.

tests/testthat/test-errors.R

@@ -223,3 +223,48 @@ test_that("unsupported complex operations are refused with R's messages", {
   }
   expect_error(quick(complex_mod), "unimplemented complex operation")
 })
+
+test_that("complex operands are refused in linear algebra", {
+  # The real BLAS/LAPACK lowerings (dgemm, dgesv, ...) would read complex
+  # storage as reals and return a plausible wrong answer where R returns a
+  # complex result: complex(2) %*% complex(2) returned a real dot product
+  # of the real parts. Refuse at compile time instead.
+  complex_matmul <- function(x, y) {
+    declare(type(x = complex(2)), type(y = complex(2)))
+    x %*% y
+  }
+  expect_error(
+    quick(complex_matmul),
+    "%*% does not support complex operands",
+    fixed = TRUE
+  )
+
+  # One complex operand is enough to poison the d* routine.
+  complex_mixed <- function(x, y) {
+    declare(type(x = complex(2, 2)), type(y = double(2, 2)))
+    x %*% y
+  }
+  expect_error(quick(complex_mixed), "does not support complex operands")
+
+  complex_solve <- function(x) {
+    declare(type(x = complex(2, 2)))
+    solve(x)
+  }
+  expect_error(quick(complex_solve), "does not support complex operands")
+
+  complex_crossprod <- function(x) {
+    declare(type(x = complex(2, 2)))
+    crossprod(x)
+  }
+  expect_error(quick(complex_crossprod), "does not support complex operands")
+
+  # t() alone is mode-preserving and keeps working on complex values.
+  complex_t <- function(x) {
+    declare(type(x = complex(2, 2)))
+    t(x)
+  }
+  expect_quick_identical(
+    complex_t,
+    list(matrix(c(1 + 1i, 2 + 0i, 3 - 1i, 4 + 2i), 2, 2))
+  )
+})