On Nov 27, 10:20�am, HardySpicer <gyansor...@gmail.com> wrote:
> On Nov 27, 4:40�am, "neophyte" <schuett@n_o_s_p_a_m.iap-kborn.de>
> wrote:
>
>
>
> > Hello,
>
> > after trying to find out why it doesn't work for several weeks, I finally
> > think I should try to find a solution with your help.
> > Here is my problem:
> > We calculate a climate modell, which consists of many array, the size of
> > the array depends on the resolution beeing choosed.
> > MAIN.f:
> > � � � implicit double precision (a-h,o-z)
> > � � � external tend
> > � � � include 'RESOLUTION' �
> > .. more variables �
> > � � � common/coplan/iplanback(maxthr),iplanforw(maxthr)
> > � � � data iplanback/maxthr*0d0/,iplanforw/maxthr*0d0/
> > c
> > c++++ destroy plans for fft routines
> > � � � call clean_plan
>
> > .. a lot of reading and pre-calculations for startup
>
> > � � � time �= 0d0
> > � � � dtime = 0.2d0*step
> > � � � do time = 0, endtime
> > � � � call tendencies(time,ym2,dydt)
> > � � � time = time + dtime
> > � � � enddo
> > ..
>
> > subroutine tend(time,y,dydt)
> > � � � USE OMP_lib
> > � � � implicit double precision (a-h,o-z)
> > � � � external intend
> > � � � include 'RESOLUTION'
> > .. other variables
>
> > !$OMP PARALLEL PRIVATE(idthr,nthr,jlat,ij,iy,w1,w2,xij)
> > !$OMP& SHARED(createbackplans,createforwplans)
> > � � � nthr=OMP_GET_NUM_THREADS()
> > � � � idthr=OMP_GET_THREAD_NUM()+1
> > !$OMP SINGLE
> > � � � numthr=nthr
> > !$OMP END SINGLE
> > !$OMP DO
> > � � � do 80 ij = 1,kh
> > � � � jlat = 2*(ij-1)+1
> > � � � xij=x(ij)
> > � � � call int(jlat,xij,w1,idthr)
> > � � � jlat = 2*(ij-1)+2
> > � � � xij=-x(ij) �
> > � � � call int(jlat,xij,w2,idthr)
> > � � � do iy = 1,ny
> > � � � work(iy,idthr) = work(iy,idthr)
> > � � �* � � � � � � � + w(ij)*(w1(iy)+w2(iy) )
> > � � � end do
> > 80 � �continue
> > !$OMP END DO
> > !$OMP END PARALLEL
> > � � � do ij=1,numthr
> > � � � do 60 iy = 1,ny
> > � 60 �dydt(iy) = dydt(iy) + work(iy,ij)
> > � � � end do
> > � � � return
> > � � � end subroutine tend
>
> > subroutine �int(jlat,x,w,idthr)
> > � � � USE OMP_lib
> > � � � implicit double precision (a-h,o-z)
> > � � � include 'RESOLUTION'
> > � � � dimension Fx(lev,0:2*nz), .more Fxs or Fys., x(irie,lev), .more xs or
> > ys.
> > .. other variables
> > � � � common/coplan/iplanback(maxthr),iplanforw(maxthr)
>
> > .. initialization of Fxs(:,:)
>
> > � � � if (iplanback(idthr) .eq. 0 ) �then
> > !$OMP Critical (createbackplans)
> > � � � � �call create_fftback_plan(Fx,idthr)
> > !$OMP End Critical (createbackplans)
> > � � � � �iplanrueck(idthr) = 1
> > � � � endif
> > � � � call fftback( Fx, x, idthr)
> > .. many more call fftback(Fx., x., idthr)
> > .. then some calculations to xs
> > � � � if (iplanforw(idthr) .eq. 0 ) �then
> > !$OMP Critical (createforwplans)
> > � � � � �call create_fftforw_plan(Fx,idthr)
> > !$OMP End Critical (createforwplans)
> > � � � � �iplanforw(idthr) = 1
> > � � � endif
> > � � � call fftforw( y, Fy, idthr)
> > .. many more call fftforw(Fy., y., idthr)
> > � � � end subroutine int
>
> > subroutine create_fftback_plan(Fx,idthr)
> > � � � USE OMP_lib
> > � � � implicit double precision (a-h,o-z)
> > � � � INCLUDE '../include/fftw3.f'
> > � � � include 'RESOLUTION'
> > � � � double precision store(lev,1:irie),ffty(lev,0:2*nz)
> > .. other variables for using fftw
> > � � � common/plancs/plan(4,maxthr)
>
> > .. some factorization to Fx and transmitting it to store
> > � � � call dfftw_plan_many_r2r(plan(3, idthr),irank,irie,multiple,
> > � � �& � � � store,inembed, istride, idist,store,onembed,
> > � � �& � � � ostride, odist, FFTW_HC2R, FFTW_MEASURE)
> > � � � end subroutine create_fftback_plan
>
> > Similar to this is the subroutine create_fftforw_plan(Fx,idthr).
>
> > subroutine fftback( Fx, x, idthr)
> > � � � USE OMP_lib
> > � � � implicit double precision (a-h,o-z)
> > � � � INCLUDE '../include/fftw3.f'
> > � � � include 'RESOLUTION'
> > � � � double precision store(lev,1:irie),ffty(lev,0:2*nz),
> > � � �& � � � � � � � � Yfftrueck(irie,lev)
> > � � � integer*8 plan
> > � � � common/plancs/plan(4,maxthr)
> > .. factorize the input and transmitting to store
> > � � � call dfftw_execute_r2r(plan(3, idthr),store,store)
> > .. factorize the output and transmitting to x
> > � � � end subroutine fftback
>
> > Again similar is subroutine fftforw( y, Fy, idthr).
>
> > And now the following problem happens. If I run the program many times with
> > only on thread, it always produces the same output. But when I start it
> > with 2 or more threads it results of the fftw-routines differ not all of
> > them but many. And because other calculations or done later, the error
> > increases. If run the program with a rough resolution, the errors are not
> > significant but with a fine resolution they do big trouble after few steps
> > I get only Zeros �or NaNs. What I figgered out already before, I was doing
> > the plannig for fftw before entering the parallel loop, but then I found
> > out for example the results from 5 threads were fine (no difference) but
> > for the other 3 threads the differed, thats why I planned for each thread
> > separately.
> > And now I have no more idea what I could or should do different.
>
> > So if anybody can help me, I would appreciate it. And sorry for the long
> > statement, I tried to put in all what could be of interest.
> > Thanks Corinna
>
> I once used Fortran - about 30 odd years ago! Bloody hell man - get
> with it!
> Climate model? I assume it must be BS in thaty case.
>
> Hardy
FORTRAN is still very commonly used for physics, and for writing C
compilers. :-)
Chris
===========================
Chris Bore
BORES SIgnal Processing
www.bores.com