external/fftpack/fftpack5/r8_factor.F

   1 subroutine r8_factor ( n, nf, fac )
   2
   3 !*****************************************************************************80
   4 !
   5 !! R8_FACTOR factors of an integer for real double precision computations.
   6 !
   7 !
   8 !
   9 !  Modified:
  10 !
  11 !    27 August 2009
  12 !
  13 !  Author:
  14 !
  15 !    Original real single precision version by Paul Swarztrauber, Dick Valent.
  16 !    Real double precision version by John Burkardt.
  17 !
  18 !  Reference:
  19 !
  20 !    Paul Swarztrauber,
  21 !    Vectorizing the Fast Fourier Transforms,
  22 !    in Parallel Computations,
  23 !    edited by G. Rodrigue,
  24 !    Academic Press, 1982.
  25 !
  26 !    Paul Swarztrauber,
  27 !    Fast Fourier Transform Algorithms for Vector Computers,
  28 !    Parallel Computing, pages 45-63, 1984.
  29 !
  30 !  Parameters:
  31 !
  32 !    Input, integer ( kind = 4 ) N, the number for which factorization and
  33 !    other information is needed.
  34 !
  35 !    Output, integer ( kind = 4 ) NF, the number of factors.
  36 !
  37 !    Output, real ( kind = 8 ) FAC(*), a list of factors of N.
  38 !
  39   implicit none
  40
  41   real ( kind = 8 ) fac(*)
  42   integer ( kind = 4 ) j
  43   integer ( kind = 4 ) n
  44   integer ( kind = 4 ) nf
  45   integer ( kind = 4 ) nl
  46   integer ( kind = 4 ) nq
  47   integer ( kind = 4 ) nr
  48   integer ( kind = 4 ) ntry
  49
  50   nl = n
  51   nf = 0
  52   j = 0
  53
  54   do while ( 1 < nl )
  55
  56     j = j + 1
  57
  58     if ( j == 1 ) then
  59       ntry = 4
  60     else if ( j == 2 ) then
  61       ntry = 2
  62     else if ( j == 3 ) then
  63       ntry = 3
  64     else if ( j == 4 ) then
  65       ntry = 5
  66     else
  67       ntry = ntry + 2
  68     end if
  69
  70     do
  71
  72       nq = nl / ntry
  73       nr = nl - ntry * nq
  74
  75       if ( nr /= 0 ) then
  76         exit
  77       end if
  78
  79       nf = nf + 1
  80       fac(nf) = real ( ntry, kind = 8 )
  81       nl = nq
  82
  83     end do
  84
  85   end do
  86
  87   return
  88 end