1      subroutine argos_pme_fft(nd1,nd2,nd3,isgn,x1,x2,
2     $     map1,map2,knode,itype)
3c $Id$
4      implicit none
5c
6#include "argos_pme_common.fh"
7#include "mafdecls.fh"
8c
9      integer nd1,nd2,nd3,isign,map1(*),map2(*),itype
10      real*8 x1(*),x2(*)
11      integer knode(ngz)
12c
13      integer ignr
14c
15c     3D-FFT interface
16c
17c     nd1,nd2,nd3  : dimension of x1 and x2
18c     isign        : 1=forward, -1=reverse, -2=reverse without normalization
19c     x1           : input slabwise distributed complex array
20c     x2           : output slabwise distributed complex array
21c     map1         : id of processor owning i-th plane
22c     map2         : id of processor owning i-th plane after transpose
23c     liwork,iwork : integer scratch space
24c     lrwork,rwork : real*8 scratch space
25c     itype        : 0=standard 3d-fft, 1=specific 3d-fft if available
26c
27      integer i,ndim,jstat,isgn
28      real*8 scale
29      integer liwrk,lrwrk,i_rwrk,l_rwrk,i_iwrk,l_iwrk
30c
31c     IBM-SP2 specific
32c
33#if defined(ESSL)
34      integer icntxt,ip(40)
35#endif
36c
37      ndim=nd1*nd2*nd3
38      isign=isgn
39      scale=one
40      if(isgn.eq.-1) scale=one/dble(ndim)
41      if(isgn.lt.-1) isign=-1
42c
43c     brand name 3d-fft routines
44c
45#if defined(ESSL) && defined(LAPI)
46      if(itype.gt.1) then
47      ip(1)=1
48      ip(2)=1
49      ip(20)=nd1
50      ip(21)=nd2
51      ip(22)=nd1
52      ip(23)=nd2
53c
54      call timer_start(206)
55      if(lpnode) call pdcft3(x1,x2,nd1,nd2,nd3,isign,scale,icntxt,ip)
56      call timer_stop(206)
57c
58      return
59      endif
60#endif
61      if(itype.gt.1)
62     + call md_abort('argos_pme_fft: 3D-pFFt no implemented',itype)
63c
64c     generic 3d-fft routine
65c
66      do 1 i=1,ngz
67      map1(i)=knode(i)
68      map2(i)=knode(i)
69    1 continue
70c
71      if(ngzloc.gt.0) then
72      do 2 i=1,nd1*nd2*mgz*2
73      x2(i)=x1(i)
74    2 continue
75c
76      call wrkspc2(1,ngx,ngy,ngz,map1,map2,liwrk,lrwrk,ignr)
77cx
78      liwrk=2*liwrk
79      lrwrk=2*lrwrk
80cx
81c
82      if(.not.ma_push_get(mt_dbl,lrwrk,'rwork',l_rwrk,i_rwrk))
83     + call md_abort('Failed to allocate rwork',0)
84c
85      if(.not.ma_push_get(mt_int,liwrk,'iwork',l_iwrk,i_iwrk))
86     + call md_abort('Failed to allocate iwork',0)
87c
88      call timer_start(206)
89      call pfft3d(isign,nd1,nd2,nd3,x2,map1,map2,
90     + lrwrk,dbl_mb(i_rwrk),liwrk,int_mb(i_iwrk),jstat)
91      call timer_stop(206)
92c
93      if(.not.ma_verify_allocator_stuff())
94     +     call md_abort('FFT buffer problems after',me)
95c      call ffflush(6)
96c
97      if(.not.ma_pop_stack(l_iwrk))
98     + call md_abort('Failed to deallocate iwork',0)
99      if(.not.ma_pop_stack(l_rwrk))
100     + call md_abort('Failed to deallocate rwork',0)
101c
102      if(isign.eq.-1) then
103      do 3 i=1,nd1*nd2*mgz*2
104      x2(i)=x2(i)*scale
105    3 continue
106      endif
107c
108      if(jstat.ne.0)
109     + call md_abort('argos_pme_fft: fft failed on node',me)
110      endif
111c
112      return
113      end
114
115