compile_sys_openmp: Compile ODE system with openmp multi-threading

In odeintr: C++ ODE Solvers Compiled on-Demand

Description

Generates an integrator using Rcpp and openmp

Usage

compile_sys_openmp(name, sys, pars = NULL, const = FALSE,
  method = "rk5_i", sys_dim = -1L, atol = 1e-06, rtol = 1e-06,
  globals = "", headers = "", footers = "", compile = TRUE,
  env = new.env(), ...)

Arguments

name	the name of the generated integration function
sys	a string containing C++ expressions
pars	a named vector of numbers or a vector of names or number of parameters
const	declare parameters const if true
method	a method string (see compile_sys)
sys_dim	length of the state vector
atol	absolute tolerance if using adaptive step size
rtol	relative tolerance if using adaptive step size
globals	a string with global C++ declarations
headers	code to appear before the odeintr namespace
footers	code to appear after the odeintr namespace
compile	if false, just return the code
env	install functions into this environment
...	passed to sourceCpp

Details

This functions behaves identically to compile_sys execpt that it does not allow one to override the default observer. In order to take advantage of openmp multi-threading, you must insert openmp pragmas into your system definition. See the examples.

A special function laplace4 is defined and can be called from your system definition. It will compute a discrete 4-point Laplacian for use in solving PDE via the method of lines. The function takes x as its first argument, dxdt as its second argument and the diffusion coefficient D as its third parameter. This function uses the default openmp scheduling.

Author(s)

Timothy H. Keitt

See Also

set_optimization, compile_sys

Examples

## Not run: 
M = 200
bistable = '
 laplace4(x, dxdt, D);  // parallel 4-point discrete laplacian
 #pragma omp parallel for
 for (int i = 0; i < N; ++i)
   dxdt[i] += a * x[i] * (1 - x[i]) * (x[i] - b);
' # bistable
compile_sys_openmp("bistable", bistable, sys_dim = M * M,
                   pars = c(D = 0.1, a = 1.0, b = 1/2),
                   const = TRUE)
at = 10 ^ (0:3)
inic = rbinom(M * M, 1, 1/2)
system.time({x = bistable_at(inic, at)})
par(mfrow = rep(2, 2), mar = rep(1, 4), oma = rep(1, 4))
for (i in 1:4){
  image(matrix(unlist(x[i, -1]), M, M),
        asp = 1, col = c("black", "lightgray"),
        axes = FALSE)
  title(main=paste("Time =", x[i, 1]))}

## End(Not run)

odeintr documentation built on May 2, 2019, 2:08 a.m.