Re: Gridding options

[Craig Markwardt <craigmnet(at)cow.physics.wisc.edu>]

tclement@ucsd.edu (Todd Clements) writes:
> arrays larger than 128x128, and starts giving incorrect values. It seems
> like just a "short" integer problem, but heck if I can figure out where it
> might be in there!!

Actually, I think this is a problem that you made NX and NY 16-bit
integers.  If you promote them to long then it should work again.

> The times are as follows (I'm not sure they mean anything for incorrect
> results, but here they are anyway):
> 
> Array size      New (1 line)       Craig (2 line)
> 512x512         0.493              0.615
> 1024x1024       2.531              3.039
> 2048x2048       10.523             12.89

Hmm, surprisingly I found that my version was about 4 times faster on
two different architectures.

          Craig      JD      1024x1024
         (otest)   (test)
Linux      0.35     1.94 s   { x86 linux unix 5.2.1 Jun  4 1999}
Alpha      1.47     6.78 s   { alpha OSF unix 5.2 Oct 30 1998}

The codes I used are below, in all their ugly, wake-up in the morning,
hair of the dog glory.  Did I do something wrong?  Note that I tried
both INT and FLOAT, and also did a comparison test.  As long as you
pass 1024L instead of 1024 you shouldn't get incorrect answers.

Craig
[cc to Clements]

.comp 
function otest, nx, iter=iter, integer=doint
ny = nx
if keyword_set(doint) then a = lindgen(nx,ny) else a = findgen(nx,ny)
if n_elements(iter) EQ 0 then iter = 10

tt = fltarr(nx+ny-1)

;; Do the work
t0 = systime(1)
ll = lindgen(nx>ny)
for j = 0, iter-1 do begin
 for i = 0, ny-1 do tt(i)      = total((a(0+ll,i-ll))(0:i<(nx-1)))
 for i = 1, nx-1 do tt(i+ny-1) = total((a(i+ll,ny-1-ll))(0:(nx-1-i)<(ny-1)))
end
if iter GT 1 then $
  print, (systime(1)-t0)/10.
return, tt
end

.comp
pro test, nx, integer=doint
ny = nx
if keyword_set(doint) then a = lindgen(nx,ny) else a = findgen(nx,ny)
t0 = systime(1)
for i = 0, 9 do $
tt=total(a[(((dy=((di=lindgen(((n=nx<ny)),nx+ny-1)))/n))*(nx gt ny?1:nx)+ $
               (nx gt ny?1:-1)*((dx=di mod n))*(nx-1))>0<(nx*ny-1)]* $
            (dy ge dx AND (dy-dx) lt nx>ny),1)
print, (systime(1)-t0)/10.
print, max(abs(tt-otest(nx,iter=1)))
end

-- 
--------------------------------------------------------------------------
Craig B. Markwardt, Ph.D.         EMAIL:    craigmnet@cow.physics.wisc.edu
Astrophysics, IDL, Finance, Derivatives | Remove "net" for better response
--------------------------------------------------------------------------