Hi every body!
This is my first day using Nim. I'm very excited by how nice and easy Nim is. After quickly skimming through the manual, I was able to translate this C code for 2D Poisson simulation: https://github.com/fortran-lang/benchmarks/blob/master/poisson2d/optimized.c into Nim.
import times
# template to simplify timed execution
template runTimed(body: untyped) =
let t1 = epochTime(); body; echo epochTime() - t1
# type
# Matrix = seq[seq[float64]]
type
Matrix[W, H: static[int]] = array[0..W-1, array[0..H-1, float64]]
proc newmat(x, y: static[int]): Matrix[x,y] =
system.zeroMem(result.addr, x*y)
# for i in 0..<x:
# for j in 0..<y:
# result[i][j] = 0.0
proc max_abs(a, b: Matrix) : float64 =
var m = 0.0
let l = a[0].len
for i in 0..<l:
for j in 0..<l:
m = max(m, abs(a[i][j]-b[i][j]))
return m
proc rho(x, y: float64) : float64 =
if 0.6 < x and x < 0.8 and 0.6 < y and y < 0.8 :
return 1.0
elif 0.2 < x and x < 0.4 and 0.2 < y and y < 0.4:
return -1.0
else :
return 0.0
proc poisson(m: static[int]) : int =
let target = 1e-6 # Target accuracy
let eps0 = 8.85e-12
let a = 0.01
let a2 = a*a
var rhoa = newmat(m,m)
var phi = newmat(m,m)
var phip = newmat(m,m)
for i in 1..m:
for j in 1..m:
rhoa[i-1][j-1] = rho(i.float64*a, j.float64*a)
# Main loop
var iter = 0
var delta = 1.0
while delta > target:
iter += 1
# Calculate new values of the potential
for i in 1..m-2:
for j in 1..m-2:
phip[i][j] = 1.0/4.0 * (phi[i+1][j] + phi[i-1][j] + phi[i][j+1] + phi[i][j-1] + a2/eps0*rhoa[i][j])
# Calculate maximum difference from old values
delta = max_abs(phip,phi)
# Swap the two arrays around
for i in 0..<m:
for j in 0..<m:
phi[i][j] = phip[i][j]
return iter
runTimed:
var c = poisson(100)
echo "m = 100, Iters = ", c
runTimed:
c = poisson(250)
echo "m = 250, Iters = ", c
Thank you.
for i in 0..<m:
for j in 0..<m:
phi[i][j] = phip[i][j]Otherwise, it looks great!
Nice code for day 1, congrats. Here's some more improvements:
type
Matrix[W, H: static[int]] = array[W, array[H, float64]]
proc initMatrix[W, H: static[int]](): Matrix[W,H] =
zeroMem(result.addr, W*H*sizeof(float64))
Matrix construction looks like initMatrix[m,m]() then.
(Note*: not completely sure about it, but arrays should be initialized to the default value of the element type automatically. So initMatrix should be obsolete and matrix construction would be just Matrix[m,m](), but that doesn't compile (why?).)
Top-level code should be encapsulated in a proc main() like this
proc main() =
runTimed:
var c = poisson(100)
echo "m = 100, Iters = ", c
runTimed:
c = poisson(250)
echo "m = 250, Iters = ", c
main()
zeromem clears a number of bytes, so the number of matrix elements has to be multiplied by the size of type float64.
Well yes it does but Nim 0 inits memory (unless otherwise told) so you can do either of the following, and it'll do the same as that zeroMem.
proc initMatrix[W, H: static[int]](): Matrix[W,H] = discard
proc initMatrix[W, H: static[int]](): Matrix[W,H] = return
Since your main solver assumes squareness (and generalizing to non-square and especially non-evenly spaced lattices changes your iterations update a lot), I might write it like this (which uses Nim features you will likely try in the near term...):
type SquareMatrix[M: static[int], F] = array[M, array[M, F]]
proc dist[M: static[int], F](a, b: SquareMatrix[M, F]): F =
for i in 0 ..< a.len: # min abs()==0==result default
for j in 0 ..< a[0].len:
result = max(result, abs(a[i][j] - b[i][j]))
proc rho[F](x, y: F) : F = # charge density
if 0.6 < x and x < 0.8 and 0.6 < y and y < 0.8: return +1.0
elif 0.2 < x and x < 0.4 and 0.2 < y and y < 0.4: return -1.0
else: return 0.0
proc poisson2Dsq[M: static[int], F](acc = 1e-6): int =
const eps0 = 8.85418781e-12 # vacuum permittivity
const a = 1.0 / float(M) # grid spacing
var a2r, phi, phip: SquareMatrix[M, F]
for i in 1..M: # populate charge density
for j in 1..M:
a2r[i-1][j-1] = a * a / eps0 * rho(i.F*a, j.F*a)
while true: # main loop
inc result # Jacobi iteration
for i in 1 ..< M-1: # update potential
for j in 1 ..< M-1:
phip[i][j] = (phi[i+1][j] + phi[i-1][j] +
phi[i][j+1] + phi[i][j-1] +a2r[i][j])/4.0
if dist(phip, phi) < acc: # calc norm(difference)
break # converged
phi = phip # save old for cvg test
You call it like poisson2Dsq[100, float32](). For me this was 3.2x faster than the float/float64 instantiation (gcc-11, nim-devel, i7-6700k), probably due to cache effects & auto-vectorization (you can pack 2x more float32's than float64's into SIMD registers & caches).
@HJarausch's advice is good if you actually need to solve these for real. This seems to come from yet another cross prog.lang benchmark suite: https://github.com/fortran-lang/benchmarks/tree/master/poisson2d
Unfortunately, it seems almost required in these to be measuring inefficient implementations/algos/ideas usually by people "new" to at least a couple of the languages in question. There is certainly a selection bias since inefficient appraoches usually amplify deltas across prog.langs, and attention is usually drawn to bigger deltas. :-(
Then use Nim's capabilities to link to an external C library.
No, use native Nim code instead. C libraries is for libs you cannot re-code yourself better in an afternoon.
I wanted to move the creation of matrices :
var phi = initMatrix[m,m]()
var phip = initMatrix[m,m]()
I'm not familiar with the problem domain but you can use var parameters for that, if needed.
proc poisson[N](a, b: var Matrix[N,N]): int =
# ...
# you can mutate a and b here
# ...
I'm very pleased to hear from you Andreas, and congrats for the great Nim. Nim community is great and very welcoming. It is my day 1 in Nim and I'm really excited by how nice, useful, and easy Nim is. Kodo for you and the Nim team.
I'm interested in scientific computing in Nim, things like Linear Algebra, Matrices, FFTs, Random Numbers, Statistcs, special functions, etc. I find Nim's syntax to be very nice, it looks like python with modern features, besides the blazing speed of C. I think it suits scientific computing tasks demanding high performance plus productivity and readability. Is it easy in Nim to interface to performance libraries (BLAS, LAPACK, FFTW, ..)? Please direct me to any similar usecases in Nim. Thank you.
No, use native Nim code instead. C libraries is for libs you cannot re-code yourself better in an afternoon.
Hmm, isn't it kinda boring to re-implement something that's already implemented? Especially if it's reasonably good and you can just use it? If the remake is a lot better or a lot faster - maybe it's worth it. But if it's just a little bit better - seems like a waste of time and effort that could be spent more productively...
Self comment:
So initMatrix should be obsolete and matrix construction would be just Matrix[m,m](), but that doesn't compile (why?).
Because it's not an object type, dummy.
Have a look at arraymancer: https://github.com/mratsim/Arraymancer
It provides wrappers to some openblas routines, and the maintainer might be interested in contributions. I added a small wrapper to *GESV a while ago (solve(...)).
More generally check out https://github.com/scinim which seems to have a wrapper for fftw3 amongst other things.
and the maintainer might be interested in contributions.
We are always happy to receive contributions or just feature requests. :)
Thanks for the mention of the SciNim organization!
In general @AboAmmar feel free to come and chat with us in the Nim Science channel, either via Matrix:
https://matrix.to/#/#nim-science:envs.net
or via the Nim Discord channel: