Merge conformability fixes into short-circuit lowering

Author: t-kalinowski

R/classes.R

@@ -339,6 +339,14 @@ Variable := new_class(
     # storage (0/1) rather than Fortran LOGICAL.
     logical_as_int = prop_bool(default = FALSE),
 
+    # Fortran kind for integer variables. User-facing R integers remain c_int;
+    # pointer-sized compiler locals opt into c_ptrdiff_t explicitly.
+    integer_kind = prop_enum(
+      c("c_int", "c_ptrdiff_t"),
+      default = "c_int",
+      exact = TRUE
+    ),
+
     # TRUE when the variable is available via host association and should not
     # be redeclared in the local scope.
     host_associated = prop_bool(default = FALSE),
@@ -350,7 +358,13 @@ Variable := new_class(
 
   validator = function(self) {
     if (isTRUE(self@logical_as_int) && !identical(self@mode, "logical")) {
-      "`logical_as_int` can only be TRUE when `mode` is 'logical'"
+      return("`logical_as_int` can only be TRUE when `mode` is 'logical'")
+    }
+    if (
+      !identical(self@integer_kind, "c_int") &&
+        !identical(self@mode, "integer")
+    ) {
+      "`integer_kind` can only differ from 'c_int' when `mode` is 'integer'"
     }
   }
 )

R/manifest.R

@@ -73,7 +73,11 @@ var_storage_bytes <- function(var) {
   switch(
     var@mode,
     double = 8,
-    integer = 4,
+    integer = if (identical(var@integer_kind, "c_ptrdiff_t")) {
+      .Machine$sizeof.pointer
+    } else {
+      4
+    },
     complex = 16,
     logical = 4,
     raw = 1,
@@ -196,7 +200,7 @@ iso_c_binding_symbols <- function(
         switch(
           var@mode,
           double = "c_double",
-          integer = "c_int",
+          integer = var@integer_kind,
           complex = "c_double_complex",
           logical = if (isTRUE(logical_is_c_int(var))) "c_int",
           raw = "c_int8_t",
@@ -262,7 +266,7 @@ emit_decl_line <- function(
   type <- switch(
     var@mode,
     double = "real(c_double)",
-    integer = "integer(c_int)",
+    integer = glue("integer({var@integer_kind})"),
     complex = "complex(c_double_complex)",
     logical = if (logical_as_int(var)) "integer(c_int)" else "logical",
     raw = "integer(c_int8_t)",
@@ -371,7 +375,7 @@ r2f.scope <- function(scope, include_errors = FALSE) {
     type <- switch(
       var@mode,
       double = "real(c_double)",
-      integer = "integer(c_int)",
+      integer = glue("integer({var@integer_kind})"),
       complex = "complex(c_double_complex)",
       logical = if (logical_as_int(var)) "integer(c_int)" else "logical",
       raw = "integer(c_int8_t)",

R/r2f-constructors.R

@@ -8,10 +8,14 @@
 # double(k), numeric(k). These lower to a single scalar literal carrying
 # array dims, so splicing contexts must spread them explicitly.
 # Used by: c(), array()
-is_fill_constructor_call <- function(e) {
-  is.call(e) &&
-    is.symbol(e[[1L]]) &&
-    as.character(e[[1L]]) %in% c("logical", "integer", "double", "numeric")
+is_fill_constructor_call <- function(e, scope) {
+  if (!is.call(e) || !is.symbol(e[[1L]])) {
+    return(FALSE)
+  }
+  name <- as.character(e[[1L]])
+  name %in%
+    c("logical", "integer", "double", "numeric") &&
+    (is.null(scope) || !inherits(scope[[name]], LocalClosure))
 }
 
 # Name of the call one frame above the current handler ("" at top level).
@@ -44,7 +48,7 @@ r2f_handlers[["c"]] <- function(args, scope = NULL, ...) {
   mode <- promoted$mode
   # Fill constructors are one scalar literal claiming length k; spread them
   # as implied-dos so the emitted element count matches the claimed length.
-  fill_idx <- which(map_lgl(args, is_fill_constructor_call))
+  fill_idx <- which(map_lgl(args, is_fill_constructor_call, scope = scope))
   if (length(fill_idx)) {
     spread_var <- NULL
     for (j in fill_idx) {
@@ -60,9 +64,16 @@ r2f_handlers[["c"]] <- function(args, scope = NULL, ...) {
           call. = FALSE
         )
       }
-      spread_var <- spread_var %||% scope_unique_var(scope, "integer")
+      spread_var <- spread_var %||%
+        scope_unique_var(
+          scope,
+          "integer",
+          integer_kind = "c_ptrdiff_t"
+        )
       ff[[j]] <- Fortran(
-        glue("({ff[[j]]}, {spread_var}=1, int({len_f}))"),
+        glue(
+          "({ff[[j]]}, {spread_var}=1_c_ptrdiff_t, int({len_f}, kind=c_ptrdiff_t))"
+        ),
         ff[[j]]@value
       )
     }
@@ -383,7 +394,7 @@ r2f_handlers[["array"]] <- function(args, scope = NULL, ..., hoist = NULL) {
       }
       shape <- glue("int([{dims_f}])")
 
-      is_fill_constructor <- is_fill_constructor_call(args$data)
+      is_fill_constructor <- is_fill_constructor_call(args$data, scope)
 
       axis_terms <- vapply(
         target_dims,

R/r2f-matrix-blas.R

@@ -61,6 +61,25 @@ assert_rhs_rank <- function(
   invisible(TRUE)
 }
 
+# BLAS/LAPACK dimensions use equality semantics: equal zero contracted
+# dimensions are conformable and can still produce a non-empty result.
+check_blas_dims <- function(left, right) {
+  if (is_wholenumber(left) && is_wholenumber(right)) {
+    return(list(
+      ok = identical(as.integer(left), as.integer(right)),
+      unknown = FALSE
+    ))
+  }
+  if (!is_scalar_na(left) && !is_scalar_na(right)) {
+    left_norm <- fortranize_expr_symbols(left)
+    right_norm <- fortranize_expr_symbols(right)
+    if (identical(left_norm, right_norm)) {
+      return(list(ok = TRUE, unknown = FALSE))
+    }
+  }
+  list(ok = TRUE, unknown = TRUE)
+}
+
 # Return the R symbol name if operand is a bare symbol; otherwise NULL.
 symbol_name_or_null <- function(x) {
   stopifnot(inherits(x, Fortran))
@@ -164,12 +183,49 @@ assert_square_matrix <- function(dims, operand, context, hoist, scope) {
     left = operand,
     right = operand,
     left_axis = 1L,
-    right_axis = 2L
+    right_axis = 2L,
+    checker = check_blas_dims
   )
 }
 
 # ---- BLAS emitters ----
 
+# Generated function results cannot currently represent zero-sized arrays.
+# Reject a known zero output during translation and guard unknown output
+# extents at runtime before emitting a BLAS call with an invalid leading
+# dimension. A zero contracted dimension remains supported when every output
+# extent is nonzero.
+assert_nonempty_blas_output <- function(
+  dim,
+  operand,
+  axis,
+  context,
+  hoist,
+  scope
+) {
+  stopifnot(
+    inherits(operand, Fortran),
+    is.numeric(axis),
+    length(axis) == 1L,
+    is_string(context)
+  )
+  message <- paste0(context, " zero-sized outputs are not supported")
+  if (is_wholenumber(dim)) {
+    if (as.integer(dim) == 0L) {
+      stop(message, call. = FALSE)
+    }
+    return(invisible(TRUE))
+  }
+
+  emit_quickr_error_if(
+    glue("{guard_dim_f(dim, operand, axis)} == 0_c_ptrdiff_t"),
+    message,
+    hoist,
+    scope
+  )
+  invisible(TRUE)
+}
+
 # Check that destination dimensions match expected output dimensions.
 assert_dest_dims_compatible <- function(dest, expected_dims, context) {
   if (is.null(dest) || is.null(expected_dims)) {
@@ -260,6 +316,32 @@ blas_int <- function(x) {
   glue("int({x_str}, kind=c_int)")
 }
 
+# Emit a BLAS call for positive contractions and fill the result with zero
+# without calling BLAS when the contracted dimension is zero.
+emit_blas_contraction <- function(call, output, contracted_dim, hoist) {
+  stopifnot(is_string(call), is_string(output))
+  assert_hoist_env(hoist)
+
+  if (is_wholenumber(contracted_dim)) {
+    if (as.integer(contracted_dim) == 0L) {
+      hoist$emit(glue("{output} = 0.0_c_double"))
+    } else {
+      hoist$emit(call)
+    }
+    return(invisible(TRUE))
+  }
+
+  hoist$emit(glue(
+    "
+if ({blas_int(contracted_dim)} == 0_c_int) then
+  {output} = 0.0_c_double
+else
+  {call}
+end if"
+  ))
+  invisible(TRUE)
+}
+
 # Centralized GEMM emission with optional destination
 # gemm: centralized BLAS GEMM emission.
 # - 'hoist' is required and provided by r2f(); handlers thread it through so
@@ -281,6 +363,22 @@ gemm <- function(
   context = "gemm"
 ) {
   assert_hoist_env(hoist)
+  assert_nonempty_blas_output(
+    m,
+    left,
+    if (opA == "N") 1L else 2L,
+    context,
+    hoist,
+    scope
+  )
+  assert_nonempty_blas_output(
+    n,
+    right,
+    if (opB == "N") 2L else 1L,
+    context,
+    hoist,
+    scope
+  )
   A_name <- ensure_blas_operand_name(left, hoist)
   B_name <- ensure_blas_operand_name(right, hoist)
 
@@ -292,18 +390,20 @@ gemm <- function(
       context = context
     )
   ) {
-    hoist$emit(glue(
+    blas_call <- glue(
       "call dgemm('{opA}','{opB}', {blas_int(m)}, {blas_int(n)}, {blas_int(k)}, 1.0_c_double, {A_name}, {blas_int(lda)}, {B_name}, {blas_int(ldb)}, 0.0_c_double, {dest@name}, {blas_int(ldc_expr)})"
-    ))
+    )
+    emit_blas_contraction(blas_call, dest@name, k, hoist)
     out <- Fortran(dest@name, dest)
     out@writes_to_dest <- TRUE
     return(out)
   }
 
   output_var <- hoist$declare_tmp(mode = "double", dims = list(m, n))
-  hoist$emit(glue(
+  blas_call <- glue(
     "call dgemm('{opA}','{opB}', {blas_int(m)}, {blas_int(n)}, {blas_int(k)}, 1.0_c_double, {A_name}, {blas_int(lda)}, {B_name}, {blas_int(ldb)}, 0.0_c_double, {output_var@name}, {blas_int(ldc_expr)})"
-  ))
+  )
+  emit_blas_contraction(blas_call, output_var@name, k, hoist)
   Fortran(output_var@name, output_var)
 }
 
@@ -325,6 +425,15 @@ gemv <- function(
   context = "gemv"
 ) {
   assert_hoist_env(hoist)
+  output_dim <- if (transA == "N") m else n
+  assert_nonempty_blas_output(
+    output_dim,
+    A,
+    if (transA == "N") 1L else 2L,
+    context,
+    hoist,
+    scope
+  )
   A_name <- ensure_blas_operand_name(A, hoist)
   x_name <- ensure_blas_operand_name(x, hoist)
 
@@ -337,18 +446,22 @@ gemv <- function(
     )
   ) {
     # Assign output to output destination
-    hoist$emit(glue(
+    blas_call <- glue(
       "call dgemv('{transA}', {blas_int(m)}, {blas_int(n)}, 1.0_c_double, {A_name}, {blas_int(lda)}, {x_name}, 1_c_int, 0.0_c_double, {dest@name}, 1_c_int)"
-    ))
+    )
+    contracted_dim <- if (transA == "N") n else m
+    emit_blas_contraction(blas_call, dest@name, contracted_dim, hoist)
     out <- Fortran(dest@name, dest)
     out@writes_to_dest <- TRUE
     return(out)
   }
   # Else assign to a temporary variable
   output_var <- hoist$declare_tmp(mode = "double", dims = out_dims)
-  hoist$emit(glue(
+  blas_call <- glue(
     "call dgemv('{transA}', {blas_int(m)}, {blas_int(n)}, 1.0_c_double, {A_name}, {blas_int(lda)}, {x_name}, 1_c_int, 0.0_c_double, {output_var@name}, 1_c_int)"
-  ))
+  )
+  contracted_dim <- if (transA == "N") n else m
+  emit_blas_contraction(blas_call, output_var@name, contracted_dim, hoist)
   Fortran(output_var@name, output_var)
 }
 
@@ -412,8 +525,6 @@ syrk <- function(
   context = "syrk"
 ) {
   assert_hoist_env(hoist)
-  X_name <- ensure_blas_operand_name(X, hoist)
-
   x_dims <- matrix_dims(X)
 
   # For trans = "T": C = t(X) %*% X, so C is k x k where k = ncol(X)
@@ -426,6 +537,15 @@ syrk <- function(
     k <- x_dims$cols
   }
   lda <- x_dims$rows
+  assert_nonempty_blas_output(
+    n,
+    X,
+    if (trans == "T") 2L else 1L,
+    context,
+    hoist,
+    scope
+  )
+  X_name <- ensure_blas_operand_name(X, hoist)
 
   # Output is symmetric n x n matrix
   writes_to_dest <- FALSE
@@ -448,9 +568,10 @@ syrk <- function(
     out_name <- out_var@name
   }
 
-  hoist$emit(glue(
+  blas_call <- glue(
     "call dsyrk('U', '{trans}', {blas_int(n)}, {blas_int(k)}, 1.0_c_double, {X_name}, {blas_int(lda)}, 0.0_c_double, {out_name}, {blas_int(n)})"
-  ))
+  )
+  emit_blas_contraction(blas_call, out_name, k, hoist)
   symmetrize_upper_to_lower(out_name, n, hoist = hoist)
 
   out <- Fortran(out_name, out_var)
@@ -547,7 +668,8 @@ triangular_solve <- function(
     left = A,
     right = B,
     left_axis = 1L,
-    right_axis = if (b_rank == 1L) NULL else 1L
+    right_axis = if (b_rank == 1L) NULL else 1L,
+    checker = check_blas_dims
   )
 
   A_name <- ensure_blas_operand_name(A, hoist)
@@ -635,7 +757,8 @@ lapack_solve <- function(
     left = A,
     right = B,
     left_axis = 1L,
-    right_axis = if (b_rank == 1L) NULL else 1L
+    right_axis = if (b_rank == 1L) NULL else 1L,
+    checker = check_blas_dims
   )
 
   A_name <- ensure_blas_operand_name(A, hoist)