Allocate block-scoped temporaries on heap

R/manifest.R

@@ -39,6 +39,63 @@ logical_as_int <- function(var) {
   identical(var@mode, "logical") && isTRUE(var@logical_as_int)
 }
 
+block_tmp_allocatable_threshold <- 16L
+
+block_tmp_element_count <- function(var) {
+  stopifnot(inherits(var, Variable))
+  dims <- var@dims
+  stopifnot(is.list(dims), length(dims) > 0L)
+  sizes <- vapply(
+    dims,
+    function(axis) {
+      if (is.integer(axis) && length(axis) == 1L && !is.na(axis)) {
+        axis
+      } else {
+        NA_integer_
+      }
+    },
+    integer(1)
+  )
+  if (anyNA(sizes)) {
+    return(NA_integer_)
+  }
+  prod(as.numeric(sizes))
+}
+
+block_tmp_allocatable <- function(
+  var,
+  scope,
+  max_stack_elements = block_tmp_allocatable_threshold
+) {
+  stopifnot(inherits(var, Variable))
+  if (!inherits(scope, "quickr_scope") || !identical(scope@kind, "block")) {
+    return(FALSE)
+  }
+  if (passes_as_scalar(var) || is.null(var@dims)) {
+    return(FALSE)
+  }
+
+  dims <- dims2f(var@dims, scope)
+  if (!nzchar(dims) || grepl(":", dims, fixed = TRUE)) {
+    return(FALSE)
+  }
+
+  n_elements <- block_tmp_element_count(var)
+  is.na(n_elements) || n_elements > max_stack_elements
+}
+
+block_tmp_allocation_lines <- function(vars, scope) {
+  stopifnot(is.list(vars))
+  allocs <- lapply(vars, function(var) {
+    if (!block_tmp_allocatable(var, scope)) {
+      return(NULL)
+    }
+    dims <- dims2f(var@dims, scope)
+    glue("allocate({var@name}({dims}))")
+  })
+  unlist(allocs, use.names = FALSE)
+}
+
 scope_vars <- function(scope) {
   vars <- as.list(scope)
   keep(vars, inherits, what = Variable)
@@ -123,15 +180,28 @@ emit_decl_line <- function(
     stop("unrecognized kind: ", format(var))
   )
 
+  # Block-scoped temporaries are explicitly marked allocatable so we can
+  # allocate them on the heap rather than relying on compiler defaults.
+  # GFortran already heap-allocates large/unknown-size locals implicitly,
+  # but flang lowers block locals to `alloca` and will stack-allocate even
+  # large runtime shapes, which can segfault under typical stack limits.
+  # We keep small, fixed-size temps (<= 16 elements) as automatic arrays
+  # to avoid allocation overhead and leave those to the compiler.
+  block_allocatable <- allow_allocatable && block_tmp_allocatable(var, scope)
+
   dims <- if (passes_as_scalar(var)) {
     NULL
+  } else if (block_allocatable) {
+    sprintf("(%s)", str_flatten_commas(rep(":", var@rank)))
   } else if (assumed_shape) {
     sprintf("(%s)", str_flatten_commas(rep(":", var@rank)))
   } else {
     dims2f(var@dims, scope) |> str_flatten_commas() |> sprintf(fmt = "(%s)")
   }
 
-  allocatable <- if (
+  allocatable <- if (block_allocatable) {
+    "allocatable"
+  } else if (
     allow_allocatable &&
       !assumed_shape &&
       !is.null(dims) &&

R/r2f.R

@@ -43,7 +43,12 @@ new_hoist <- function(scope) {
     stmts <- str_split_lines(hoisted, code)
 
     if (has_block()) {
-      decls <- emit_decls(scope_vars(block_scope), block_scope)
+      block_vars <- scope_vars(block_scope)
+      decls <- emit_decls(block_vars, block_scope)
+      allocs <- block_tmp_allocation_lines(block_vars, block_scope)
+      if (length(allocs)) {
+        stmts <- c(allocs, stmts)
+      }
       return(str_flatten_lines(emit_block(decls, stmts)))
     }
 

doc/flang-block-temporaries.md

@@ -0,0 +1,154 @@
+# Flang, `block`-scoped temporaries, and stack vs heap allocation
+
+This note documents why quickr uses Fortran `block` scopes for temporaries, what we observed when compiling quickr-generated Fortran with different compilers (flang vs gfortran), and the practical trade-offs for users and maintainers.
+
+## Why quickr uses `block` for temporaries
+
+quickr emits intermediate temporaries (e.g. work buffers and copies of inputs for LAPACK calls) into a Fortran `block` so their lifetime is limited:
+
+- It reduces *peak* live memory within a subroutine by constraining temporary lifetimes to the smallest region that needs them.
+- It gives the compiler clearer lifetime boundaries (enter/exit) so storage can be reclaimed as early as possible.
+
+Example shape (simplified):
+
+```fortran
+block
+  real(c_double) :: tmp(m, n)
+  ! ... use tmp ...
+end block
+! tmp is out of scope here
+```
+
+## Minimal example and generated Fortran
+
+This minimal quickr function forces quickr to create runtime-sized temporaries in a `block` (rectangular least-squares solve):
+
+```r
+f_rect <- function(X, y) {
+  declare(type(X = double(n, k)), type(y = double(n)))
+  solve(X, y)
+}
+```
+
+The generated Fortran (see `scratch/solve_rect.f90`) includes automatic arrays with runtime bounds inside a `block`, e.g.:
+
+```fortran
+block
+  real(c_double) :: btmp1_(X__dim_1_, X__dim_2_)
+  real(c_double) :: btmp2_(max(X__dim_1_, X__dim_2_), 1)
+  ! ...
+end block
+```
+
+These are *explicit local variables* (not compiler-invented temporaries), and their bounds are only known at runtime.
+
+## What flang does: stack allocation via `alloca`
+
+When we compile the generated Fortran with flang to LLVM IR:
+
+```sh
+flang-new -S -emit-llvm -O0 scratch/solve_rect.f90 -o scratch/solve_rect.ll
+```
+
+the IR shows `block` entry/exit implemented as a stack save/restore, with runtime-sized locals allocated on the stack:
+
+```llvm
+%34 = call ptr @llvm.stacksave.p0()
+%45 = alloca double, i64 %44, align 8
+%52 = alloca double, i64 %51, align 8
+%64 = alloca double, i64 %63, align 8
+call void @llvm.stackrestore.p0(ptr %34)
+```
+
+Interpretation:
+
+- The `block` lifetime is respected (stack space is released at `stackrestore`).
+- But the storage comes from the thread stack. Large `m*n` can exceed OS stack limits and segfault.
+
+The flang assembly (`scratch/solve_rect_flang.s`) also reflects stack allocation: it computes an aligned size and adjusts `sp` (stack pointer) down by that amount, rather than calling `malloc`.
+
+## What gfortran does: heap allocation via `malloc/free`
+
+When we compile the same Fortran with gfortran:
+
+```sh
+gfortran -S -O0 scratch/solve_rect.f90 -o scratch/solve_rect_gfortran.s
+```
+
+the assembly shows heap allocation and cleanup for these runtime-sized locals:
+
+- Calls to `_malloc` appear at the points where runtime-sized arrays are created.
+- Calls to `_free` appear at the end of the scope/cleanup path.
+
+Example snippet shape (see `scratch/solve_rect_gfortran.s`):
+
+```asm
+bl _malloc
+...
+bl _malloc
+...
+bl _malloc
+...
+bl _free
+bl _free
+bl _free
+```
+
+This behavior matches the common gfortran strategy of heap-allocating large/unknown-size locals, which avoids stack overflows at the cost of allocator overhead.
+
+## Why flang flags didn’t help here
+
+`flang-new --help` exposes:
+
+- `-fno-stack-arrays` (“Allocate array temporaries on the heap (default)”)
+- `-fstack-arrays` (“Attempt to allocate array temporaries on the stack, no matter their size”)
+
+However, these apply to *compiler-generated array temporaries*. In this case, quickr emits explicit local automatic arrays in the source:
+
+```fortran
+real(c_double) :: btmp1_(X__dim_1_, X__dim_2_)
+```
+
+Recompiling `scratch/solve_rect.f90` with `-fno-stack-arrays` still produced dynamic `alloca` in LLVM IR for these locals, i.e. flang continued to allocate them on the stack.
+
+## Trade-offs and implications
+
+**Stack allocation (flang behavior here)**
+
+- Pros: fast allocation, fast reclamation (especially with `block`), low fragmentation.
+- Cons: brittle for large runtime sizes; depends on OS thread stack limits; may crash abruptly (segfault) rather than produce a clean R error.
+
+**Heap allocation (gfortran behavior here)**
+
+- Pros: can handle much larger work arrays without stack overflow; less sensitive to stack limits.
+- Cons: allocator overhead; potential fragmentation; requires reliable cleanup on all exit paths.
+
+**`block` is still useful**
+
+Even when heap-allocating, `block` remains a good design for minimizing peak live memory: it constrains lifetimes so temporaries can be freed early (by the compiler/runtime).
+
+## User guidance: choosing a different Fortran compiler
+
+quickr selects a Fortran compiler based on a single global option:
+
+```r
+options(quickr.fortran_compiler = "gfortran")  # force gfortran
+options(quickr.fortran_compiler = "flang")     # force flang
+options(quickr.fortran_compiler = "auto")      # default behavior
+```
+
+Notes:
+
+- On macOS, `auto` currently prefers flang if it is available.
+- If flang compilation fails, quickr falls back to gfortran and disables automatic flang preference for the rest of the session.
+- When OpenMP is used, quickr may avoid flang unless it was explicitly requested.
+
+## Recommended direction for quickr (future work)
+
+If we want `block`-scoped temporaries *and* robust behavior across compilers:
+
+- Avoid emitting large runtime-sized temporaries as automatic arrays.
+- Emit them as `allocatable` locals within the `block` and `allocate(...)` them explicitly, so they live on the heap while still having short lifetimes.
+
+This keeps the original motivation for `block` (reduced peak memory) while avoiding flang’s stack allocation behavior for large runtime-sized locals.
+

tests/testthat/_snaps/block-scopes.md

@@ -77,3 +77,89 @@
         return out;
       }
 
+# block-scoped temps allocate on the heap for runtime shapes
+
+    Code
+      cat("# Snapshot note: ", note, "\n", sep = "")
+    Output
+      # Snapshot note: Block temps with runtime sizes are allocatable so flang doesn't stack-allocate large work arrays.
+    Code
+      fn
+    Output
+      function(x) {
+          declare(type(x = double(n, m)))
+          out <- ifelse((x > 0.0)[1, 1], 1.0, 0.0)
+          out
+        }
+      <environment: 0x0>
+    Code
+      cat(fsub)
+    Output
+      subroutine fn(x, out, x__dim_1_, x__dim_2_) bind(c)
+        use iso_c_binding, only: c_double, c_int
+        implicit none
+
+        ! manifest start
+        ! sizes
+        integer(c_int), intent(in), value :: x__dim_1_
+        integer(c_int), intent(in), value :: x__dim_2_
+
+        ! args
+        real(c_double), intent(in) :: x(x__dim_1_, x__dim_2_)
+        real(c_double), intent(out) :: out
+        ! manifest end
+
+
+        block
+          logical, allocatable :: btmp1_(:, :) ! logical
+
+          allocate(btmp1_(x__dim_1_, x__dim_2_))
+          btmp1_ = ((x > 0.0_c_double))
+          out = merge(1.0_c_double, 0.0_c_double, btmp1_(1_c_int, 1_c_int))
+        end block
+      end subroutine
+    Code
+      cat(cwrapper)
+    Output
+      #define R_NO_REMAP
+      #include <R.h>
+      #include <Rinternals.h>
+
+
+      extern void fn(
+        const double* const x__, 
+        double* const out__, 
+        const R_len_t x__dim_1_, 
+        const R_len_t x__dim_2_);
+
+      SEXP fn_(SEXP _args) {
+        // x
+        _args = CDR(_args);
+        SEXP x = CAR(_args);
+        if (TYPEOF(x) != REALSXP) {
+          Rf_error("typeof(x) must be 'double', not '%s'", Rf_type2char(TYPEOF(x)));
+        }
+        const double* const x__ = REAL(x);
+        const int* const x__dim_ = ({
+        SEXP dim_ = Rf_getAttrib(x, R_DimSymbol);
+        if (Rf_length(dim_) != 2) Rf_error(
+          "x must be a 2D-array, but length(dim(x)) is %i",
+          (int) Rf_length(dim_));
+        INTEGER(dim_);});
+        const int x__dim_1_ = x__dim_[0];
+        const int x__dim_2_ = x__dim_[1];
+
+        const R_xlen_t out__len_ = (1);
+        SEXP out = PROTECT(Rf_allocVector(REALSXP, out__len_));
+        double* out__ = REAL(out);
+
+        fn(
+          x__,
+          out__,
+          x__dim_1_,
+          x__dim_2_);
+
+        UNPROTECT(1);
+        return out;
+      }
+

tests/testthat/_snaps/closure-hoist-snapshots.md

@@ -49,8 +49,9 @@
 
 
             block
-              real(c_double) :: btmp1_(x__len_)
+              real(c_double), allocatable :: btmp1_(:)
 
+              allocate(btmp1_(x__len_))
               btmp1_ = ((x + 1.0_c_double))
               res = btmp1_(i)
             end block
@@ -139,14 +140,16 @@
 
 
             block
-              real(c_double) :: btmp1_(nx, ny)
+              real(c_double), allocatable :: btmp1_(:, :)
 
+              allocate(btmp1_(nx, ny))
               btmp1_ = ((temp + 1.0_c_double))
               temp(1_c_int, 1_c_int) = btmp1_(1_c_int, 1_c_int)
             end block
             block
-              real(c_double) :: btmp1_(nx, ny)
+              real(c_double), allocatable :: btmp1_(:, :)
 
+              allocate(btmp1_(nx, ny))
               btmp1_ = ((temp + 2.0_c_double))
               temp(nx, ny) = btmp1_(nx, ny)
             end block