Good morning,
I am struggling to speed up my fortran code. The problem seems simple but I don't know if it is really possible to speed up these loops in a simple way without using a specific toolbox or parallelization and I didn't find an answer for know (I hope it is not too obvious otherwise i am sorry...).
Here my problem:
I would like to convert a scattered distribution of dense points [xs,ys,zs,rs] on a coarser regular grid taking the mean of the values inside each cell defined by vectors [xv,yv,zv]. This requires several loop to test each grid block and it takes a very long time to run even for relatively small grid size and I have to scale up at least to 10 times larger grid in each direction compared to this example.
Do somebody have an idea how to speed up this kind of code ? or any hint on how I should proceed to improve my code.
Thank you. Gamma
EDIT: if soemone has another method to pass from a scattered point distribution to a regular grid it is also nice thanks
Here's a code sample:
PROGRAM build_scatter2reg_gridmodel
IMPLICIT NONE
REAL, DIMENSION(:), ALLOCATABLE :: xs,ys,zs,rs
REAL, DIMENSION(:), ALLOCATABLE :: rsloc
REAL :: rsloc2
REAL, DIMENSION(:), ALLOCATABLE :: xv,yv,zv
REAL, DIMENSION(:,:,:), ALLOCATABLE :: xg,yg,zg,rg
INTEGER :: i,j,k,l
INTEGER :: nx,ny,nz
INTEGER, DIMENSION(:), ALLOCATABLE :: idex
REAL, DIMENSION(3) :: xmax,xmin
REAL :: dx,dy,dz
CHARACTER(LEN=100) :: filename
INTEGER :: lu,answ
INTEGER :: irec,nlines
INTEGER :: N
character(8) :: date
character(10) :: time
character(5) :: zone
integer,dimension(8) :: values
lu=66
nlines=300000
!allocate vector
ALLOCATE(xs(nlines),ys(nlines),zs(nlines),rs(nlines))
ALLOCATE(rsloc(nlines))
xmin(1)=-1000
xmin(2)=-5000
xmin(3)=0
xmax(1)=1000
xmax(2)=5000
xmax(3)=5000
!random fill
call random_number(xs)
xs = xmin(1) + (xmax(1)-xmin(1))*xs
call random_number(ys)
ys = xmin(2) + (xmax(2)-xmin(2))*ys
call random_number(zs)
zs = xmin(3) + (xmax(3)-xmin(3))*zs
call random_number(rs)
rs = 0 + (10-0)*rs
nx=20
ny=20
nz=20
dx= (xmax(1) - xmin(1))/(nx-1)
dy= (xmax(2) - xmin(2))/(ny-1)
dz= (xmax(3) - xmin(3))/(nz-1)
call date_and_time(date,time)
WRITE(*,*) 'date start:', date,' , ', time
ALLOCATE(xv(nx),yv(ny),zv(nz))
xv(:) = (/ ( xmin(1) + (i-1)*dx ,i=1,nx)/)
yv(:) = (/ ( xmin(2) + (i-1)*dy ,i=1,ny)/)
zv(:) = (/ ( xmin(3) + (i-1)*dz ,i=1,nz)/)
!OPEN (lu, file ='output_build_scatter.asc')
open(unit=lu,file='output_build_scatter.asc',form='unformatted',access='direct', &
status='replace',recl=4)
!Fill and write 3d grid
irec=0
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
!using COUNT/WHERE
N=COUNT( xs .ge. xv(i) .AND. xs .lt. xv(i+1) .AND. &
ys .ge. yv(j) .AND. ys .lt. yv(j+1) .AND. &
zs .ge. zv(k) .AND. zs .lt. zv(k+1) )
rsloc=0
WHERE( xs .ge. xv(i) .AND. xs .lt. xv(i+1) .AND. &
ys .ge. yv(j) .AND. ys .lt. yv(j+1) .AND. &
zs .ge. zv(k) .AND. zs .lt. zv(k+1) )
rsloc= rs/real(N)
!ELSEWHERE
ENDWHERE
WRITE(lu,rec=1+irec) xv(i)
WRITE(lu,rec=2+irec) yv(i)
WRITE(lu,rec=3+irec) zv(i)
WRITE(lu,rec=4+irec) sum(rsloc)
irec=irec+4
ENDDO
ENDDO
WRITE(*,*) 'i,j,k',i,j,k
ENDDO
CLOSE (lu)
call date_and_time(date,time)
WRITE(*,*) 'date end:', date,' , ', time
END PROGRAM