Add private variables to openMP loops

[mns-nordicals]

Added support for using private variable in openMP loops. The implementation adds the new api:

• declare(parallel(private = x))) - for single private variables
• declare(parallel(private = c(x, y))) - for multiple private variables.

The private variables must be in scope or we error. This is to protect against putting a wrong variable name in private.

I tried to add code that tracks if a scalar is changed and error if it is not private. This is in order to not create race conditions and possibly compile errors. However, this might be a bit difficult and there is likely many edge cases that cannot be anticipated. So I will leave it up to you how much the user should be protected from wrong use of openMP.

"Motivating example"

roll_mean <- function(x, window) {
  declare(
    type(x = double(n)),
    type(window = integer(1))
  )

  out <- double(length(x) - window + 1)

  declare(parallel())
  for (i in seq_along(out)) {
    acc <- 0
    for (j in seq_len(window)) {
      acc <- acc + x[(i + j - 1)]
    }
    out[i] <- acc / window
  }
  out
}

This is a simplified version of the roll_mean example. An LLM suggested this would be faster than just slapping the openMP declaration over the original example - turns out it didn't make much difference. But at least it is possible now.

[mns-nordicals]

@codex review

[chatgpt-codex-connector]

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: a4b381e50a

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

• Open a pull request for review
• Mark a draft as ready
• Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

[chatgpt-codex-connector]

Use scoped Fortran names for OpenMP private list

When a private variable is shadowed in a closure scope, its Fortran symbol can be renamed via make_shadow_fortran_name (e.g., i__local_). Here the private list is built with fortranize_name, which ignores those scope-specific renames and emits the base name instead. In nested closures or when a variable name collides with a parent scope, this produces !$omp parallel do private(i) even though the actual loop variable is i__local_, leading to a compile error (“symbol not declared”) or leaving the true variable non-private. Consider deriving private names from the scope (e.g., scope_fortran_symbol or get0(name, scope)@name) rather than re-fortranizing the raw symbol.

Useful? React with 👍 / 👎.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 92.81%. Comparing base (f22ae80) to head (4f72349).
⚠️ Report is 57 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #88      +/-   ##
==========================================
+ Coverage   92.76%   92.81%   +0.05%     
==========================================
  Files          27       27              
  Lines        5734     5777      +43     
==========================================
+ Hits         5319     5362      +43     
  Misses        415      415              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[t-kalinowski]

Thanks for putting this together!

I was under the (perhaps mistaken!) impression that OpenMP and most compilers will automatically infer block locals as thread private variables, if the block is declared within an OpenMP region. Is that not the case?

[mns-nordicals]

Current situation on main

So currently my example on main generates this code:

[note: I have reference example at the bottom]

Loop version

subroutine roll_mean_loop(x, window, out, x__len_) bind(c)
  use iso_c_binding, only: c_double, c_int, c_ptrdiff_t
  implicit none

  ! manifest start
  ! sizes
  integer(c_ptrdiff_t), intent(in), value :: x__len_

  ! args
  real(c_double), intent(in) :: x(x__len_)
  integer(c_int), intent(in) :: window
  real(c_double), intent(out) :: out(((x__len_ - window) + 1))

  ! locals
  integer(c_int) :: i
  real(c_double) :: acc
  integer(c_int) :: j
  ! manifest end


  out = 0

  !$omp parallel do
  do i = 1, size(out)
    acc = 0.0_c_double
    do j = 1, window
      acc = (acc + x((((i + j) - 1_c_int))))
    end do
    out(i) = (acc / real(window, kind=c_double))
  end do
  !$omp end parallel do
end subroutine

As I understand it, the initial loop index i is automatically local by openMP, however that , j acc are global. This can create race conditions where two threads modify the same variable.

sapply version

While an sapply version creates a local subroutine that gets called

subroutine roll_mean_sapply(x, window, out, x__len_) bind(c)
  use iso_c_binding, only: c_double, c_int, c_ptrdiff_t
  implicit none

  ! manifest start
  ! sizes
  integer(c_ptrdiff_t), intent(in), value :: x__len_

  ! args
  real(c_double), intent(in) :: x(x__len_)
  integer(c_int), intent(in) :: window
  real(c_double), intent(out) :: out(((x__len_ - window) + 1))

  ! locals
  integer(c_int) :: out_len
  integer(c_int) :: tmp1_
  ! manifest end



  out_len = ((size(x) - window) + 1_c_int)
  out = 0

  !$omp parallel do
  do tmp1_ = 1_c_int, ((x__len_ - window) + 1_c_int)
    call roll_mean(tmp1_, out(tmp1_))

  end do
  !$omp end parallel do

  contains
    subroutine roll_mean(i, res)
      use iso_c_binding, only: c_double, c_int
      implicit none

      integer(c_int), intent(in) :: i
      real(c_double), intent(out) :: res
      real(c_double) :: acc
      integer(c_int) :: j

      acc = 0.0_c_double
      do j = 1, window
        acc = (acc + x((((i + j) - 1_c_int))))
      end do
      res = (acc / real(window, kind=c_double))
    end subroutine
end subroutine

So as I understand it, the sapply version is safe as i, j acc and res are local to the roll_mean subroutine.

Summary of situation on main

The loop version does not create a local block like:

!$omp parallel do
do i = 1, size(out)
  block
    real(c_double) :: acc
    integer(c_int) :: j
    acc = 0.0_c_double
    do j = 1, window
      acc = acc + x(i + j - 1)
    end do
    out(i) = acc / real(window, kind=c_double)
  end block
end do
!$omp end parallel do

This is also an alternative route to my PR. I think one of the difficulties is to infer what should be local and what should be global. For instance out should not be local.

I will note that, when I run my two example below, that I get the same results, even though they are supposedly modifying the same variable. So maybe the compiler is somehow inferring the locals.

NOTE ! after consulting an AI it seems my example is somewhat poorly constructed Claude/chatGPT had this to say about the safety of my example on main

• acc - It's assigned (acc = 0.0_c_double) before being read in each iteration of the outer loop. Many compilers (including gfortran) can detect this pattern.
• j - As the loop variable of the inner do loop, gfortran treats it as private because it's a loop index within the parallel region.

and

Most importantly for your concern:

• acc must be private for correctness.
  • In OpenMP, variables that are assigned within the loop body and are not meant to be shared typically need to be private.
  • For a parallel do, compilers implement the OpenMP standard’s data-sharing rules: a variable like acc that is used as a temporary inside the loop must behave as thread-local (otherwise it would be a blatant race and would break a huge amount of real code).

So the correct results you’re seeing are because each thread has its own acc and j (and its own i), not because the compiler “inferred” something ad-hoc.

My PR

With declare(parallel(private = c(acc, j))) it generates this code on the loop example

subroutine roll_mean_loop(x, window, out, x__len_) bind(c)
  use iso_c_binding, only: c_double, c_int, c_ptrdiff_t
  implicit none

  ! manifest start
  ! sizes
  integer(c_ptrdiff_t), intent(in), value :: x__len_

  ! args
  real(c_double), intent(in) :: x(x__len_)
  integer(c_int), intent(in) :: window
  real(c_double), intent(out) :: out(((x__len_ - window) + 1))

  ! locals
  integer(c_int) :: i
  real(c_double) :: acc
  integer(c_int) :: j
  ! manifest end


  out = 0

  !$omp parallel do private(acc, j)
  do i = 1, size(out)
    acc = 0.0_c_double
    do j = 1, window
      acc = (acc + x((((i + j) - 1_c_int))))
    end do
    out(i) = (acc / real(window, kind=c_double))
  end do
  !$omp end parallel do
end subroutine

I was actually not fully aware how the sapply version worked, so i think the last commit can be roll back - my sapply example at least now fails with the private declaration ...

Summary

I think the sapply version is probably fine - but I think we need to change something with a loop openMP version. Unless I am specifying something wrong for the lock scope not to be created.

Reference example

Loop openMP

roll_mean_loop <- function(x, window) {
  declare(
    type(x = double(n)),
    type(window = integer(1))
  )

  out <- double(length(x) - window + 1)

  declare(parallel()) # and my PR has private = c(acc, j) inside parallel()
  for (i in seq_along(out)) {
    acc <- 0
    for (j in seq_len(window)) {
      acc <- acc + x[(i + j - 1)]
    }
    out[i] <- acc / window
  }
  out
}

sapply openMP

roll_mean_sapply <- function(x, window) {
  declare(
    type(x = double(n)),
    type(window = integer(1))
  )

  roll_mean <- function(i) {
    acc <- 0
    for (j in seq_len(window)) {
      acc <- acc + x[(i + j - 1)]
    }
    acc / window
  }

  out_len <- length(x) - window + 1L
  out <- double(out_len)

  declare(parallel())
  out <- sapply(seq_len(out_len), roll_mean)

  out
}

[mns-nordicals]

I apologies for all these long comments.

The loop version of openMP still needs some work.

Actually unsafe example

I tried to create a new example that is actually unsafe - a parallel sum.

library(quickr)

loop_sum_unsafe <- function(x) {
  declare(type(x = double(n)))

  total <- 0.0

  declare(parallel())
  for (i in seq_along(x)) {
    total <- total + x[i]  # Race! All threads read-modify-write `total`
  }

  total
}

loop_sum_safe <- function(x) {
  declare(type(x = double(n)))

  total <- 0.0

  declare(parallel(private = total))
  for (i in seq_along(x)) {
    total <- total + x[i]  
  }

  total
}

x <- rnorm(1000, 5, 5)

unsafe_parallel_sum <- quick(loop_sum_unsafe)
safe_parallel_sum <- quick(loop_sum_safe)

quickr:::r2f(loop_sum_unsafe)
quickr:::r2f(loop_sum_safe)

r_loop_sum <- loop_sum_unsafe

r_loop_sum(x)
unsafe_parallel_sum(x)
safe_parallel_sum(x)

When running this (ignoring the generated code) I get

[1] 5028.085
[1] 631.6976
[1] 0

So at least this demonstrates unsafe openMP loop without private declaration.

It also show my implementation is lacking as I now declared total as private so each thread get their own version and what is return is the global one that I set to 0.

There seem to be many openMP clauses that you can set. In this case, it seem like we need reduction(+:total)

So maybe we should just close this and think how to extend openmp loops - i didn't test with the sapply version, maybe it works. Maybe sapply version should be the only version.

**** ADDED
After reading a bit more about openMP I see that it is quite more flexible and with more options than I first realized.

Do you want to limit quickr's openMP capabilities to map/apply style functions where input size and output size matches? In that case as noted above, I think allowing openMP with for-loop might be too flexible a prone to errors.

[mns-nordicals]

Needs more thought - but it would be nice to support more OpenMP declarations.

[t-kalinowski]

Apologies for not responding yet! This is a very thoughtful issue that identifies a real need, and I’m still collecting my thoughts on how to respond. It’s still on my to-do list to write up an issue describing how I’m imagining this interface will evolve. We will need a way to expose private variables and other syntax that can go into an OpenMP directive, but before doing so, I want to lay the groundwork for how declare(parallel()) might work with OpenACC and CUDA as well, and then figure out the common ground between those and OpenMP-specific directives.

[t-kalinowski]

I wouldn’t normally share raw output from ChatGPT, but in the interest of time, here is the last turn from a brainstorming ChatGPT session I had a few weeks ago, while thinking about how to evolve the parallel interface:

ChatGPT output

Based on the tile/block-first patterns described in the cuTile summary you shared (hide threads, expose tiles, give smart defaults, make “kernel-ish” code feel like normal array code), here’s an updated, concrete declare() design for CPU parallelism + GPU offload in quickr.

I’m aiming for:

• Naive users: “make this loop run in parallel / on GPU” (one word).
• Power users: pass through OpenMP/OpenACC/CUDA-ish clauses (minimal new syntax).
• A cuTile-like mental model: users specify tile/block size and let the compiler handle the rest.

---

1) Core attachment model

declare() stores compile-time annotations that apply to the next eligible construct, then are consumed.

Eligible constructs (initially):

• for (i in <index-iterable>) { ... }
• out <- sapply(seq_len(n), FUN) and out <- sapply(seq_along(x), FUN) (your existing lowering)

Also add one more target:

• a following { ... } block (for data/target regions)

---

2) Minimal surface API

A. The one-liner for naive users

declare(parallel())
for (i in seq_len(n)) { ... }

declare(gpu())
for (i in seq_len(n)) { ... }

• parallel() means “CPU parallelize next loop”.
• gpu() means “offload next loop to GPU” (backend chosen by compiler/config).

Both are “friendly defaults”.

B. A cuTile-inspired knob: tile/block size

declare(gpu(tile = 256))
for (i in seq_len(n)) { ... }

Interpretation: “run on GPU, chunk the iteration space into tiles of ~256 iterations per block”.

• For CUDA-ish backends this maps naturally to threads per block.
• For OpenACC/OpenMP offload this maps to a sensible vector length / thread_limit hint.

Also allow 2D tiling for collapsed loops:

declare(gpu(tile = c(16, 16), collapse = 2))
for (i in seq_len(n)) {
  for (j in seq_len(m)) {
    ...
  }
}

C. Data/target regions (avoid repeated copies)

A block-scoped data region attaches to the next { ... } block:

declare(gpu_data(copyin = vars(a, b), copyout = vars(y)))
{
  declare(gpu(tile = 256))
  for (i in seq_len(n)) y[i] <- a[i] + b[i]

  declare(gpu(tile = 256))
  for (i in seq_len(n)) y[i] <- y[i] * 2
}

• This is the closest quickr analogue of cuTile’s “global arrays live on device for the kernel”.
• You still annotate loops with gpu() (explicit is better than surprising), but the data mapping is shared.

D. “vars()” helper for symbol lists

vars(a, b, y)

This stays unevaluated inside declare() and is easy to validate.

---

3) Unifying CPU and GPU: parallel() as the generic entry point

If you want one concept that scales from CPU to GPU:

declare(parallel(device = "cpu"))
declare(parallel(device = "gpu"))

With a backend override:

declare(parallel(device = "gpu", backend = "acc"))  # OpenACC
declare(parallel(device = "gpu", backend = "omp"))  # OpenMP target offload
declare(parallel(device = "gpu", backend = "cuda")) # CUDA-ish lowering

And keep the shorthands:

• gpu(...) := parallel(device="gpu", ...)
• cpu(...) := parallel(device="cpu", ...) (optional)
• omp(...), acc(...), cuda(...) as explicit backend shorthands (below)

---

4) Backend-specific shorthands for power users

These are still passed via declare(), but allow more direct clause control.

OpenMP (CPU or target offload)

declare(omp())
for (i in seq_len(n)) { ... }

Common knobs:

declare(omp(
  device      = "cpu",      # or "gpu" for target offload
  num_threads = 8,
  schedule    = "static",
  collapse    = 2,
  private     = vars(tmp)
))
for (...) { ... }

OpenACC (GPU)

declare(acc())
for (i in seq_len(n)) { ... }

With tiling-ish intent:

declare(acc(tile = 256))     # maps to vector_length / vector
for (...) { ... }

CUDA-ish directive over a loop

declare(cuda(kernel = TRUE, tile = 256))
for (i in seq_len(n)) { ... }

This is the “cuTile vibe”: you say “kernel + tile size”, and don’t talk about threads explicitly.

---

5) Reductions in an R-like way

Make this easy and safe:

declare(gpu(reduction = sum(total)))
for (i in seq_len(n)) {
  total <- total + x[i]
}

Multiple reductions:

declare(gpu(reduction = list(sum = total, max = best)))
for (i in seq_len(n)) {
  total <- total + x[i]
  best  <- max(best, x[i])
}

Backend mapping examples:

• OpenMP: reduction(+:total) / reduction(max:best)
• OpenACC: reduction(+:total) etc
• CUDA-ish: lower to warp/block reductions if supported; otherwise reject or fall back

(You can start by supporting sum() only, then expand.)

---

6) Pass-through escape hatch

This keeps experts unblocked.

declare(acc(raw = "parallel loop gang vector_length(256)"))
for (i in seq_len(n)) { ... }

declare(omp(raw = "target teams distribute parallel do collapse(2)"))
for (...) { ... }

declare(cuda(raw = "cuf kernel do(1) <<<*,256>>>"))
for (...) { ... }

If a raw directive is not a structured form, optionally allow:

declare(acc(raw = "...", end = "..."))

---

7) Defaulting rules inspired by cuTile

These defaults make the “one word” API actually work:

declare(gpu()) defaults

• tile: "auto" (compiler chooses; maybe 256 for 1D loops)
• data: "auto" for single-loop offload
  (compiler infers read-only vs write arrays and emits per-loop mapping)
• backend: "auto" (pick from available toolchain; user can override via option/env var)

declare(gpu_data(...)) defaults

• If a var is in copyout, it’s assumed written on device.
• If omitted, gpu_data() with no args could mean “present/managed/auto” (depending on backend choice).

---

8) Concrete examples

A) Naive GPU vector add

f <- quick(function(a, b) {
  declare(
    type(a = double(n)),
    type(b = double(n)),
    type(y = double(n))
  )

  declare(gpu())
  for (i in seq_len(n)) {
    y[i] <- a[i] + b[i]
  }

  y
})

B) cuTile-like “tile size” tuning

declare(gpu(tile = 256))
for (i in seq_len(n)) {
  y[i] <- a[i] + b[i]
}

C) Keep data on GPU across loops

declare(gpu_data(copyin = vars(a, b), copyout = vars(y)))
{
  declare(gpu(tile = 256))
  for (i in seq_len(n)) y[i] <- a[i] + b[i]

  declare(gpu(tile = 256))
  for (i in seq_len(n)) y[i] <- y[i] * 2
}

D) sapply() lowering to a GPU loop

declare(gpu())
out <- sapply(seq_len(n), \(i) x[i] * 2)

(Your existing compile_sapply_assignment() becomes the insertion point.)

---

9) What I would ship first

If you want a staged rollout that still “feels cuTile-ish” immediately:

1. declare(gpu()) on a for loop
   Emits: OpenACC parallel loop (or CUDA cuf kernel do) with auto data mapping.
2. tile = <int> as a single tuning knob.
3. gpu_data(...) applied to {...} blocks (explicit regions).
4. Reductions: reduction = sum(x) (only) to start.
5. Backend shorthands + raw= pass-through.

That gives you a simple story early (“GPU loop + tile size + auto data”), while leaving room for the full OpenACC/OpenMP/CUDA clause richness later.

If you want, I can also translate this into a precise grammar (what calls/args are accepted as language objects) and a mapping table showing how each knob expands to OpenACC vs OpenMP target vs CUDA-ish directives in generated Fortran.

---

About 70% of that strikes me as the wrong direction, but hopefully it gives a sense of where I’m heading with this.

Before evolving toward GPU support, quickr needs a design/interface for transferring data to and from the GPU. We don’t want to force a data move every time we enter and exit a parallel region, which means quickr probably needs to evolve toward being able to define “kernel” functions that take a pointer to data that already lives on the GPU.

That means the next step is adding the ability to accept non-atomic types as inputs and return non-atomic types as outputs. A reasonable initial goal might be to support defining a quickr function that takes a NumPy array (via reticulate). That would involve exposing an interface for other R packages to register code generators for their “type” (e.g., an S3 class for numpy.ndarray). These would need to generate R code, C-bridge code, and Fortran manifest and signature code. Also, this would require thinking through what the declare(type(...)) syntax would be for these external pointer types. This is not a trivial task, but I’d like to tackle it before expanding the parallel interface.

On a slightly orthogonal note, I think quickr already has enough information to automatically infer which variables are private and/or what the reduction operation is. We should be able to generate the appropriate OpenMP directive automatically, without requiring the user to specify it.

[t-kalinowski]

It would be good to fix the correctness bug you identified however. I suspect that forcing parallel loop bodies to always emig BLOCK is the easiest approach.

[mns-nordicals]

Thank you for your reply and for providing details on your future direction. For what you want to support I see that this requires careful thinking and planning.

For what it's worth, I think an api that looks something like this:

declare(parallel(
  num_threads = 8,
  schedule    = "static",
  reduction = sum(total), 
  private     = vars(tmp)
))
for (...) { ... }

looks interesting and aligns with my own thoughts.

This gives a lot of flexibility but off-loads a lot of responsibility to the user. If it is possible to auto infer some stuff like reductions and private variables that would be nice and create opportunity of a user writing declare(parallel()) for a "just works" experience.