program hpftest3
c
c This program calls the MPL library spfft. 
c Example using advanced initialization.
c
       integer nxbits, nybits
       parameter (nxbits=7,nybits=7)
c  128 by 128 MasPar
       integer dimx, dimy
       parameter (dimx=512,dimy=1024)
c  NB: dimx and dimy must be multiples of nxproc and nyproc, respectively
       complex a(0:dimx-1, 0:dimy-1)
       integer ctrl
       integer i, trszx , trszy, mb(0:2), agrlst(0:3)

C prepare for doing 32x32 transforms stored "hierarchical" on 128x128 maspar
C 9 and 10 is the size of the array (2^9 x 2^10)
C There must always be a bit in mb(2)
       mb(0) = 9-nxbits
       mb(1) = 10-nybits-1
       mb(2) = 1
       nactgr = 2
       agrlst(0) = 0
c size of transform 2^5 in x-dir
       agrlst(1) = 5
c 9 is the first y-direction bit
       agrlst(2) = 9
c size of transform 2^5 in y-dir
       agrlst(3) = 5

c MPL subroutines
CMPF MPL  fft_ainit, sp_fft, fft_free test

c Initialize
       forall (i=0:dimx-1,j=0:dimy-1) a(i,j) = cmplx(i,-j)


c      call fft_ainit (ctrl,mode,mb,nactgr,agrlst,postprmut)
c          in the following call we do not have a post-permutation
       call fft_ainit (ctrl, 0, mb, nactgr, agrlst, 0) 

       write(*,*) 'Input:'
       write(*,*) 'a(0,1) =', a(0,1)
       write(*,*) 'a(97,17) =', a(97,17)
       write(*,*) 'a(200,250) =', a(200,250)

c      Forward 2-D transform
       write(*,*) 'After FFT:'
       call sp_fft(a,1,ctrl)
       write(*,*) 'a(0,1) =', a(0,1)
       write(*,*) 'a(97,17) =', a(97,17)
       write(*,*) 'a(200,250) =', a(200,250)


c      Inverse 2-D transform
       call sp_fft(a,-1,ctrl)

       write(*,*) 'After inverse FFT:'
       write(*,*) 'a(0,1) =', a(0,1)
       write(*,*) 'a(97,17) =', a(97,17)
       write(*,*) 'a(200,250) =', a(200,250)

c      Free temporary data structures
       call fft_free(ctrl)

       end