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