Re: Do I need MathLink to run finite-difference fast enough for

*To*: mathgroup at smc.vnet.net*Subject*: [mg115887] Re: Do I need MathLink to run finite-difference fast enough for*From*: Oliver Ruebenkoenig <ruebenko at wolfram.com>*Date*: Mon, 24 Jan 2011 05:23:13 -0500 (EST)

On Sun, 23 Jan 2011, Oliver Ruebenkoenig wrote: > On Sat, 22 Jan 2011, James wrote: > >> Hello Mike, >> >> This is a Brusselator model of chemical reaction diffussion that is >> employing the Euler numerical method to solve the PDE system: >> >> u_t=D_u lap(u)+a-(b+1)u+u^2v >> >> v_t=D_v lap(v)+bu-u^2v >> >> For suitable values of the parameters, the characteristic Turing >> patterns of dots and stripes should emerge. >> >> Using Oliver's Compiler suggestions above under V7, I can run the >> simulator on a 64X64 grid 10000 times on my machine in about 5.5 >> seconds. >> However, I feel that's still too slow to create a reasonable interactive >> Demonstration of the Brusselator where the user will be changing the >> paramters. Ideally, I'd like to get it down to about 2 seconds. I'm >> not >> sure how fast it would run in V8. Here is the complete code using >> Oliver's suggestions with some additional improvements for stripes that >> runs in 5.5 seconds on my machine. Can anyone suggest of a way to get >> it down to two seconds? >> >> n = 64; >> a = 4.5; >> b = 7.5; >> du = 2; >> dv = 16; >> dt = 0.01; >> totaliter = 10000; >> u = a + 0.3*RandomReal[{-0.5, 0.5}, {n, n}]; >> v = b/a + 0.3*RandomReal[{-0.5, 0.5}, {n, n}]; >> cf = Compile[{{uIn, _Real, 2}, {vIn, _Real, 2}, >> {aIn, _Real}, {bIn, _Real}, {duIn, _Real}, >> {dvIn, _Real}, {dtIn, _Real}, {iterationsIn, >> _Integer}}, Block[{u = uIn, v = vIn, lap, dt = dtIn, >> k = bIn + 1, kern = N[{{0, 1, 0}, {1, -4, 1}, >> {0, 1, 0}}], du = duIn, dv = dvIn}, >> Do[lap = RotateLeft[u, {1, 0}] + RotateLeft[u, >> {0, 1}] + RotateRight[u, {1, 0}] + >> RotateRight[u, {0, 1}] - 4*u; >> u = u + dt*(du*lap + a - u*(k - v*u)); >> lap = RotateLeft[v, {1, 0}] + RotateLeft[v, >> {0, 1}] + RotateRight[v, {1, 0}] + >> RotateRight[v, {0, 1}] - 4*v; >> v = v + dt*(dv*lap + u*(b - v*u)); , >> {iterationsIn}]; {u, v}]]; >> Timing[c1 = cf[u, v, a, b, du, dv, dt, totaliter]; ] >> ListDensityPlot[c1[[1]]] >> >> > > There are two omissions in the above code. a should be aIn and b should > be bIn. > > OK, this fixed code runs on my laptop in V7 in about 9.8 sec. > > In V8 this runs in about 6.5 sec. If I then use CompilationTarget -> "C" > it runs in about 2.9 sec. (actually in a second run it then takes about > 2.6 sec, because the M- and the compiler then have initialized everything > they need) > > > > > James, here is a further small speedup: cf = Compile[{{uIn, _Real, 2}, {vIn, _Real, 2}, {aIn, _Real}, {bIn, _Real}, {duIn, _Real}, {dvIn, _Real}, \ {dtIn, _Real}, {iterationsIn, _Integer}}, Block[{u = uIn, v = vIn, lap, dt = dtIn, k = bIn + 1, kern = N@{{0, 1, 0}, {1, -4, 1}, {0, 1, 0}}, du = duIn, dv = dvIn}, Do[ lap = RotateLeft[u, {1, 0}] + RotateLeft[u, {0, 1}] + RotateRight[u, {1, 0}] + RotateRight[u, {0, 1}] - 4*u; u += dt*(du*lap + aIn - u*(k - v*u)); lap = RotateLeft[v, {1, 0}] + RotateLeft[v, {0, 1}] + RotateRight[v, {1, 0}] + RotateRight[v, {0, 1}] - 4*v; v += dt*(dv*lap + u*(bIn - v*u)); , {iterationsIn}]; {u, v}], CompilationTarget -> "C", RuntimeOptions -> "Speed" ]; this changes u = u + ... to u+=... and sets RuntimeOptions->"Speed". I also tried to parallelize this, but the approach I choose had little success: cflap = Compile[{{uIn, _Real, 1}}, RotateLeft[uIn] + RotateRight[uIn], RuntimeAttributes -> {Listable}, CompilationTarget -> "C", RuntimeOptions -> "Speed",Parallelization -> True] re = cflap[u] - 4*u + Transpose[cflap[Transpose[u]]]; I guess here the overhead is too big. Oliver

**Follow-Ups**:**Re: Vector problem***From:*Mark Perrin <m.perrin@me.com>

**Vector problem***From:*Mark Perrin <m.perrin@me.com>