| How to speed up KRIG2D by 30x [message #86116] |
Tue, 08 October 2013 07:51  |
![Mike[5] is currently offline](theme/default/images/xoffline.png.pagespeed.ic.XRkd1fkXye.png "Mike[5] is currently offline")
Mike[5]
Messages: 5
Registered: October 2013
|
Junior Member |
|
|
While trying to figure out what the KRIG2D routine was doing, I discovered it can be speed up by one order of magnitude or more, for large numbers of output points.
One only needs to change a single line of code.
I post this here hoping someone at EXELIS will see it and include in the next IDL distribution...
In practice the last few lines of the current KRIG2D routine (as of IDL 8.2.3) should be changed from this
;++++++++++++++++++++++++++++++
for j=0,ny-1 do begin ;Each output point
y0 = bounds[1] + gs[1] * j
for i=0,nx-1 do begin
x0 = bounds[0] + gs[0] * i
d[0] = reform(sqrt((x-x0)^2 + (y-y0)^2),n) ;distance
d = call_function(fname, d, t) ;Get rhs
d[n] = 1.0 ;lagrange constr
d=LUSOL( a, indx, d) ;Solve using LUDC results.
r[i,j] = total(d * z)
endfor
endfor
;------------------------------
into the following mathematically equivalent form, with LUSOL outside the double FOR loop and a change in the order of the array operations. This works because array multiplication is commutative.
;++++++++++++++++++++++++++++++
az = LUSOL(a, indx, [z,0.0]) ;Solve using LUDC results.
for j=0,ny-1 do begin ;Each output point
y0 = bounds[1] + gs[1] * j
for i=0,nx-1 do begin
x0 = bounds[0] + gs[0] * i
d[0] = reform(sqrt((x-x0)^2 + (y-y0)^2),n) ;distance
d = call_function(fname, d, t) ;Get rhs
d[n] = 1.0 ;lagrange constr
r[i,j] = total(d * az)
endfor
endfor
;------------------------------
In the following example, with a large number of output points, my modified routine KRIG2D_MODIFIED is 30x faster than IDL's KRIG2D and also 7x faster than the native code GRIDDATA. This suggests that GRIDDATA is making the same mistake as KRIG2D of including the back substitution inside the final for loop.
;++++++++++++++++++++++++++++++
; Create a dataset of N points.
n = 500 ;# of scattered points
seed = -121147L ;For consistency
x = RANDOMU(seed, n)
y = RANDOMU(seed, n)
; Create a dependent variable in the form a function of (x,y)
z = 3 * EXP(-((9*x-2)^2 + (7-9*y)^2)/4) + $
3 * EXP(-((9*x+1)^2)/49 - (1-0.9*y)) + $
2 * EXP(-((9*x-7)^2 + (6-9*y)^2)/4) - $
EXP(-(9*x-4)^2 - (2-9*y)^2)
s = [0.5,0]
tic
z2 = krig2d_modified(z, x, y, EXPONENTIAL=s, NX=200, NY=200)
toc
tic
z2 = GRIDDATA(x, y, z, /KRIGING, DIMENSION=[200,200], VARIOGRAM=[2,s])
toc
tic
z2 = krig2d(z, x, y, EXPONENTIAL=s, NX=200, NY=200)
toc
;------------------------------
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
