program test1
C
      parameter (NBMAX = 20, NMAX = 2**NBMAX)
      real x(0:NMAX-1)
      integer rhs (NBMAX), mat (NBMAX * NBMAX)
      integer nbits, i, err
      integer A, R

C     MPL subroutines in p3pack library
CMPF MPL blp_btprmNew, blp_rout, pmr_blpCheck, permut32, blp_print

C     Use MPF directives to ensure the control arrays are stored on front end.
C     Note that this is currently the default unless you use the arrays in 
C     HPF syntax. Thus, in this example the directive could have been omitted.
CMPF ONFE rhs, mat
C     Store the actual data array on the DPU.
CMPF ONDPU x

      write (*,10) 'Size of problem (try e.g. 16, max is ', NBMAX, '): '
      read (*,*) nbits
  10  format (A,I2,A)
      write (*,*) 'Trying a perfect shuffle on an array of length', 
     &            2**nbits

C Make perfect shuffle permutation
      do i = 1, nbits
         rhs (i) = i
      enddo
      rhs (nbits) = 0
      
C Set up bit permutation
      call blp_btprmNew (A, rhs, nbits)

C Print the permutation
      write (*,*) 'Here is the permutation:'
      call blp_print (A)

C Compute routing information
      call blp_rout (R, A)

C Check if everything worked OK
      call pmr_blpCheck (err, A, R)

      if (err.gt.0) then
         write (*, *) 'Error in permutation!'
      else
         write (*, *) 'Permutation is correct!'
      endif

C Initialize data
      forall (i= 0:2**nbits-1) x(i) = real(i)

      write(*,*) 
      write(*,*) 'Input:'
      write(*,*) 'x(0) =', x(0)
      write(*,*) 'x(1) =', x(1)
      write(*,*) 'x(456) =', x(456)
      write(*,*) 'x(457) =', x(457)
      write(*,*) 'x(max) =', x(2**nbits-1)

C Do the actual permutation on the data. We can use either permut or permut32.
C permut32 is slightly faster.
C     call permut(x, 4, 1, 1, R)
      call permut32(x, 1, R)

C Check answer
      write(*,*) 
      write(*,*) 'Output: (after perfect shuffle)'
      write(*,*) 'x(0) =', x(0)
      write(*,*) 'x(1) =', x(1)
      write(*,*) 'x(456) =', x(456)
      write(*,*) 'x(457) =', x(457)
      write(*,*) 'x(max) =', x(2**nbits-1)

      end