******************************************************************
  
  subroutine pda_cfftf(n,c,wsave)
  
  ******************************************************************
  
  subroutine pda_cfftf computes the forward complex discrete fourier
  transform (the fourier analysis). equivalently , pda_cfftf computes
  the fourier coefficients of a complex periodic sequence.
  the transform is defined below at output parameter c.
  
  the transform is not normalized. to obtain a normalized transform
  the output must be divided by n. otherwise a call of pda_cfftf
  followed by a call of pda_cfftb will multiply the sequence by n.
  
  the array wsave which is used by subroutine pda_cfftf must be
  initialized by calling subroutine pda_cffti(n,wsave).
  
  input parameters
  
  
  n      the length of the complex sequence c. the method is
         more efficient when n is the product of small primes. n
  
  c      a complex array of length n which contains the sequence
  
  wsave   a real work array which must be dimensioned at least 4n+15
          in the program that calls pda_cfftf. the wsave array must be
          initialized by calling subroutine pda_cffti(n,wsave) and a
          different wsave array must be used for each different
          value of n. this initialization does not have to be
          repeated so long as n remains unchanged thus subsequent
          transforms can be obtained faster than the first.
          the same wsave array can be used by pda_cfftf and pda_cfftb.
  
  output parameters
  
  c      for j=1,...,n
  
             c(j)=the sum from k=1,...,n of
  
                   c(k)*exp(-i*(j-1)*(k-1)*2*pi/n)
  
                         where i=sqrt(-1)
  
  wsave   contains initialization calculations which must not be
          destroyed between calls of subroutine pda_cfftf or pda_cfftb