Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Ignore:
Timestamp:
Jul 6, 2020 3:56:08 PM (4 years ago)
Author:
raasch
Message:

files re-formatted to follow the PALM coding standard

File:
1 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/singleton_mod.f90

    r4182 r4591  
    11!> @file singleton_mod.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
     3! This file is part of the PALM model system.
     4!
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
     15!
     16! Copyright 1997-2020 Leibniz Universitaet Hannover
     17!--------------------------------------------------------------------------------------------------!
     18!
     19!
    320! Current revisions:
    421! -----------------
    5 ! 
    6 ! 
     22!
     23!
    724! Former revisions:
    825! -----------------
    926! $Id$
     27! File re-formatted to follow the PALM coding standard
     28!
     29!
     30! 4182 2019-08-22 15:20:23Z scharf
    1031! Corrected "Former revisions" section
    11 ! 
     32!
    1233! 3761 2019-02-25 15:31:42Z raasch
    13 ! statement added to prevent compiler warning about unused variables
     34! Statement added to prevent compiler warning about unused variables
    1435!
    1536! Revision 1.1  2002/05/02 18:56:59  raasch
    … …  
    1738!
    1839!
     40!--------------------------------------------------------------------------------------------------!
    1941! Description:
    2042! ------------
    2143!> Multivariate Fast Fourier Transform
    2244!>
    23 !> Fortran 90 Implementation of Singleton's mixed-radix algorithm,
    24 !> RC Singleton, Stanford Research Institute, Sept. 1968.
    25 !>
    26 !> Adapted from fftn.c, translated from Fortran 66 to C by Mark Olesen and
    27 !> John Beale.
    28 !>
    29 !> Fourier transforms can be computed either in place, using assumed size
    30 !> arguments, or by generic function, using assumed shape arguments.
    31 !>
     45!> Fortran 90 Implementation of Singleton's mixed-radix algorithm, RC Singleton, Stanford Research
     46!> Institute, Sept. 1968.
     47!>
     48!> Adapted from fftn.c, translated from Fortran 66 to C by Mark Olesen and John Beale.
     49!>
     50!> Fourier transforms can be computed either in place, using assumed size arguments, or by generic
     51!> function, using assumed shape arguments.
    3252!>
    3353!> Public:
    3454!>
    35 !>   fftkind                              kind parameter of complex arguments
    36 !>                                        and function results.
     55!>   fftkind                              kind parameter of complex arguments and function results.
    3756!>
    3857!>   fft(array, dim, inv, stat)           generic transform function
    … …  
    5271!> Formal Parameters:
    5372!>
    54 !>   array    The complex array to be transformed. array can be of arbitrary
    55 !>            rank (i.e. up to seven).
    56 !>
    57 !>   shape    With subroutine fftn, the shape of the array to be transformed
    58 !>            has to be passed separately, since fftradix - the internal trans-
    59 !>            formation routine - will treat array always as one dimensional.
    60 !>            The product of elements in shape must be the number of
     73!>   array    The complex array to be transformed. Array can be of arbitrary rank (i.e. up to seven).
     74!>
     75!>   shape    With subroutine fftn, the shape of the array to be transformed has to be passed
     76!>            separately, since fftradix - the internal transformation routine - will always treat
     77!>            array as one dimensional. The product of elements in shape must be the number of
    6178!>            elements in array.
    62 !>            Although passing array with assumed shape would have been nicer,
    63 !>            I prefered assumed size in order to prevent the compiler from
    64 !>            using a copy-in-copy-out mechanism. That would generally be
    65 !>            necessary with fftn passing array to fftradix and with fftn
    66 !>            being prepared for accepting non consecutive array sections.
    67 !>            Using assumed size, it's up to the user to pass an array argu-
    68 !>            ment, that can be addressed as continous one dimensional array
    69 !>            without copying. Otherwise, transformation will not really be
    70 !>            performed in place.
    71 !>            On the other hand, since the rank of array and the size of
    72 !>            shape needn't match, fftn is appropriate for handling more than
    73 !>            seven dimensions.
    74 !>            As far as function fft is concerned all this doesn't matter,
    75 !>            because the argument will be copied anyway. Thus no extra
    76 !>            shape argument is needed for fft.
     79!>            Although passing array with assumed shape would have been nicer, I prefered assumed
     80!>            size in order to prevent the compiler from using a copy-in-copy-out mechanism. That
     81!>            would generally be necessary with fftn passing array to fftradix and with fftn being
     82!>            prepared for accepting non consecutive array sections. Using assumed size, it's up to
     83!>            the user to pass an array argument, that can be addressed as continous one dimensional
     84!>            array without copying. Otherwise, transformation will not really be performed in place.
     85!>            On the other hand, since the rank of array and the size of shape needn't match, fftn
     86!>            is appropriate for handling more than seven dimensions. As far as function fft is
     87!>            concerned all this doesn't matter, because the argument will be copied anyway. Thus no
     88!>            extra shape argument is needed for fft.
    7789!>
    7890!> Optional Parameters:
    7991!>
    80 !>   dim      One dimensional integer array, containing the dimensions to be
    81 !>            transformed. Default is (/1,...,N/) with N being the rank of
    82 !>            array, i.e. complete transform. dim can restrict transformation
    83 !>            to a subset of available dimensions. Its size must not exceed the
    84 !>            rank of array or the size of shape respectivly.
    85 !>
    86 !>   inv      If .true., inverse transformation will be performed. Default is
    87 !>            .false., i.e. forward transformation.
    88 !>
    89 !>   stat     If present, a system dependent nonzero status value will be
    90 !>            returned in stat, if allocation of temporary storage failed.
     92!>   dim      One dimensional integer array, containing the dimensions to be transformed. Default
     93!>            is (/1,...,N/) with N being the rank of array, i.e. complete transform. dim can
     94!>            restrict transformation to a subset of available dimensions. Its size must not exceed
     95!>            the rank of array or the size of shape respectivly.
     96!>
     97!>   inv      If .true., inverse transformation will be performed. Default is .false., i.e. forward
     98!>            transformation.
     99!>
     100!>   stat     If present, a system dependent nonzero status value will be returned in stat, if
     101!>            allocation of temporary storage failed.
    91102!>
    92103!>
    93104!> Scaling:
    94105!>
    95 !>   Transformation results will always be scaled by the square root of the
    96 !>   product of sizes of each dimension in dim. (See examples below)
     106!>   Transformation results will always be scaled by the square root of the product of sizes of each
     107!>   dimension in dim. (See examples below)
    97108!>
    98109!>
    … …  
    111122!>     result = fft(A, dim=(/1,3/))
    112123!>
    113 !>   will transform with respect to the first and the third dimension, scaled
    114 !>   by sqrt(L*N).
     124!>   will transform with respect to the first and the third dimension, scaled by sqrt(L*N).
    115125!>
    116126!>     result = fft(fft(A), inv=.true.)
    … …  
    127137!>
    128138!>   Following changes have been introduced with respect to fftn.c:
    129 !>   - complex arguments and results are of type complex, rather than
    130 !>     real an imaginary part separately.
    131 !>   - increment parameter (magnitude of isign) has been dropped,
    132 !>     inc is always one, direction of transform is given by inv.     
    133 !>   - maxf and maxp have been dropped. The amount of temporary storage
    134 !>     needed is determined by the fftradix routine. Both fftn and fft
    135 !>     can handle any size of array. (Maybe they take a lot of time and
    136 !>     memory, but they will do it)
    137 !>
    138 !>   Redesigning fftradix in a way, that it handles assumed shape arrays
    139 !>   would have been desirable. However, I found it rather hard to do this
    140 !>   in an efficient way. Problems were:
    141 !>   - to prevent stride multiplications when indexing arrays. At least our
    142 !>     compiler was not clever enough to discover that in fact additions
    143 !>     would do the job as well. On the other hand, I haven't been clever
    144 !>     enough to find an implementation using array operations.
    145 !>   - fftradix is rather large and different versions would be necessaray
    146 !>     for each possible rank of array.
    147 !>   Consequently, in place transformation still needs the argument stored
    148 !>   in a consecutive bunch of memory and can't be performed on array
    149 !>   sections like A(100:199:-3, 50:1020). Calling fftn with such sections
    150 !>   will most probably imply copy-in-copy-out. However, the function fft
    151 !>   works with everything it gets and should be convenient to use.
     139!>   - Complex arguments and results are of type complex, rather than real an imaginary part
     140!>     separately.
     141!>   - Increment parameter (magnitude of isign) has been dropped, inc is always one, direction of
     142!>     transform is given by inv.
     143!>   - maxf and maxp have been dropped. The amount of temporary storage needed is determined by the
     144!>     fftradix routine. Both fftn and fft can handle any size of array. (Maybe they take a lot of
     145!>     time and memory, but they will do it)
     146!>
     147!>   Redesigning fftradix in a way, that it handles assumed shape arrays would have been desirable.
     148!>   However, I found it rather hard to do this in an efficient way. Problems were:
     149!>   - To prevent stride multiplications when indexing arrays. At least our compiler was not clever
     150!>     enough to discover that in fact additions would do the job as well. On the other hand, I
     151!>     haven't been clever enough to find an implementation using array operations.
     152!>   - fftradix is rather large and different versions would be necessaray for each possible rank of
     153!>     array.
     154!>   Consequently, in place transformation still needs the argument stored in a consecutive bunch of
     155!>   memory and can't be performed on array sections like A(100:199:-3, 50:1020). Calling fftn with
     156!>   such sections will most probably imply copy-in-copy-out. However, the function fft works with
     157!>   everything it gets and should be convenient to use.
    152158!>
    153159!> Michael Steffens, 09.12.96, <Michael.Steffens@mbox.muk.uni-hannover.de>
    154160!> Restructured fftradix for better optimization. M. Steffens, 4 June 1997
    155 !------------------------------------------------------------------------------!
     161!--------------------------------------------------------------------------------------------------!
    156162 MODULE singleton
    157  
     163
    158164
    159165    USE kinds
    … …  
    162168
    163169    PRIVATE
    164     PUBLIC:: fft, fftn
    165 
    166     REAL(wp), PARAMETER:: sin60 = 0.86602540378443865_wp
    167     REAL(wp), PARAMETER:: cos72 = 0.30901699437494742_wp
    168     REAL(wp), PARAMETER:: sin72 = 0.95105651629515357_wp
    169     REAL(wp), PARAMETER:: pi    = 3.14159265358979323_wp
     170    PUBLIC ::  fft   !<
     171    PUBLIC ::  fftn  !<
     172
     173    REAL(wp), PARAMETER ::  cos72 = 0.30901699437494742_wp  !<
     174    REAL(wp), PARAMETER ::  pi    = 3.14159265358979323_wp  !<
     175    REAL(wp), PARAMETER ::  sin60 = 0.86602540378443865_wp  !<
     176    REAL(wp), PARAMETER ::  sin72 = 0.95105651629515357_wp  !<
    170177
    171178    INTERFACE fft
    … …  
    183190
    184191
    185 !------------------------------------------------------------------------------!
     192!--------------------------------------------------------------------------------------------------!
    186193! Description:
    187194! ------------
    188195!> @todo Missing function description.
    189 !------------------------------------------------------------------------------!
    190  FUNCTION fft1d(array, dim, inv, stat) RESULT(ft)
     196!--------------------------------------------------------------------------------------------------!
     197 FUNCTION fft1d( array, dim, inv, stat ) RESULT( ft )
    191198!
    192199!-- Formal parameters
    193     COMPLEX(wp),  DIMENSION(:), INTENT(IN)           :: array
    194     INTEGER(iwp), DIMENSION(:), INTENT(IN),  OPTIONAL:: dim
    195     INTEGER(iwp),               INTENT(OUT), OPTIONAL:: stat
    196     LOGICAL,                    INTENT(IN),  OPTIONAL:: inv
     200    COMPLEX(wp),  DIMENSION(:), INTENT(IN)            ::  array  !<
     201    INTEGER(iwp), DIMENSION(:), INTENT(IN),  OPTIONAL ::  dim    !<
     202    INTEGER(iwp),               INTENT(OUT), OPTIONAL ::  stat   !<
     203    LOGICAL,                    INTENT(IN),  OPTIONAL ::  inv    !<
    197204!
    198205!-- Function result
    199     COMPLEX(wp), DIMENSION(SIZE(array, 1)):: ft
    200 
    201     INTEGER(iwp)::  idum
    202     INTEGER(iwp)::  ishape(1)
     206    COMPLEX(wp), DIMENSION(SIZE(array, 1)) ::  ft  !<
     207
     208    INTEGER(iwp) ::  idum       !<
     209    INTEGER(iwp) ::  ishape(1)  !<
    203210
    204211!
    … …  
    208215    ft = array
    209216    ishape = SHAPE( array )
    210     CALL fftn(ft, ishape, inv = inv,  stat = stat)
     217    CALL fftn( ft, ishape, inv = inv,  stat = stat )
    211218!
    212219!-- Next statement to prevent compiler warning about unused variable
    … …  
    216223
    217224
    218 !------------------------------------------------------------------------------!
     225!--------------------------------------------------------------------------------------------------!
    219226! Description:
    220227! ------------
    221228!> @todo Missing function description.
    222 !------------------------------------------------------------------------------!
    223  FUNCTION fft2d(array, dim, inv, stat) RESULT(ft)
     229!--------------------------------------------------------------------------------------------------!
     230 FUNCTION fft2d( array, dim, inv, stat ) RESULT( ft )
    224231!
    225232!-- Formal parameters
    226     COMPLEX(wp),  DIMENSION(:,:), INTENT(IN)           :: array
    227     INTEGER(iwp), DIMENSION(:),   INTENT(IN),  OPTIONAL:: dim
    228     INTEGER(iwp),                 INTENT(OUT), OPTIONAL:: stat
    229     LOGICAL,                      INTENT(IN),  OPTIONAL:: inv
     233    COMPLEX(wp),  DIMENSION(:,:), INTENT(IN)            ::  array  !<
     234    INTEGER(iwp), DIMENSION(:),   INTENT(IN),  OPTIONAL ::  dim    !<
     235    INTEGER(iwp),                 INTENT(OUT), OPTIONAL ::  stat   !<
     236    LOGICAL,                      INTENT(IN),  OPTIONAL ::  inv    !<
    230237!
    231238!-- Function result
    232     COMPLEX(wp), DIMENSION(SIZE(array, 1), SIZE(array, 2)):: ft
    233 
    234     INTEGER(iwp) ::  ishape(2)
     239    COMPLEX(wp), DIMENSION(SIZE(array, 1), SIZE(array, 2)) ::  ft  !<
     240
     241    INTEGER(iwp) ::  ishape(2)  !<
     242!
     243!-- Intrinsics used
     244    INTRINSIC SIZE, SHAPE
     245
     246    ft = array
     247    ishape = SHAPE( array )
     248    CALL fftn( ft, ishape, dim, inv, stat )
     249
     250 END FUNCTION fft2d
     251
     252
     253!--------------------------------------------------------------------------------------------------!
     254! Description:
     255! ------------
     256!> @todo Missing function description.
     257!--------------------------------------------------------------------------------------------------!
     258 FUNCTION fft3d( array, dim, inv, stat ) RESULT( ft )
     259!
     260!-- Formal parameters
     261    COMPLEX(wp),  DIMENSION(:,:,:), INTENT(IN)            ::  array  !<
     262    INTEGER(iwp), DIMENSION(:),     INTENT(IN),  OPTIONAL ::  dim    !<
     263    INTEGER(iwp),                   INTENT(OUT), OPTIONAL ::  stat   !<
     264    LOGICAL,                        INTENT(IN),  OPTIONAL ::  inv    !<
     265!
     266!-- Function result
     267    COMPLEX(wp), DIMENSION(SIZE(array, 1), SIZE(array, 2), SIZE(array, 3)) ::  ft  !<
     268
     269    INTEGER(iwp) ::  ishape(3)  !<
     270
     271!
     272!-- Intrinsics used
     273    INTRINSIC SIZE, SHAPE
     274
     275    ft = array
     276    ishape = SHAPE( array)
     277    CALL fftn(ft, ishape, dim, inv, stat)
     278
     279 END FUNCTION fft3d
     280
     281
     282!--------------------------------------------------------------------------------------------------!
     283! Description:
     284! ------------
     285!> @todo Missing function description.
     286!--------------------------------------------------------------------------------------------------!
     287 FUNCTION fft4d( array, dim, inv, stat ) RESULT( ft )
     288!
     289!-- Formal parameters
     290    COMPLEX(wp),  DIMENSION(:,:,:,:), INTENT(IN)            ::  array  !<
     291    INTEGER(iwp), DIMENSION(:),       INTENT(IN),  OPTIONAL ::  dim    !<
     292    INTEGER(iwp),                     INTENT(OUT), OPTIONAL ::  stat   !<
     293    LOGICAL,                          INTENT(IN),  OPTIONAL ::  inv    !<
     294!
     295!-- Function result
     296    COMPLEX(wp), DIMENSION(SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4)) ::  ft  !<
     297
     298    INTEGER(iwp) ::  ishape(4)  !<
    235299!
    236300!-- Intrinsics used
    … …  
    241305    CALL fftn(ft, ishape, dim, inv, stat)
    242306
    243  END FUNCTION fft2d
    244 
    245 
    246 !------------------------------------------------------------------------------!
     307 END FUNCTION fft4d
     308
     309
     310!--------------------------------------------------------------------------------------------------!
    247311! Description:
    248312! ------------
    249313!> @todo Missing function description.
    250 !------------------------------------------------------------------------------!
    251  FUNCTION fft3d(array, dim, inv, stat) RESULT(ft)
     314!--------------------------------------------------------------------------------------------------!
     315 FUNCTION fft5d( array, dim, inv, stat ) RESULT( ft )
    252316!
    253317!-- Formal parameters
    254     COMPLEX(wp),  DIMENSION(:,:,:), INTENT(IN)           :: array
    255     INTEGER(iwp), DIMENSION(:),     INTENT(IN),  OPTIONAL:: dim
    256     INTEGER(iwp),                   INTENT(OUT), OPTIONAL:: stat
    257     LOGICAL,                        INTENT(IN),  OPTIONAL:: inv
     318    COMPLEX(wp),  DIMENSION(:,:,:,:,:), INTENT(IN)            ::  array  !<
     319    INTEGER(iwp), DIMENSION(:),         INTENT(IN),  OPTIONAL ::  dim    !<
     320    INTEGER(iwp),                       INTENT(OUT), OPTIONAL ::  stat   !<
     321    LOGICAL,                            INTENT(IN),  OPTIONAL ::  inv    !<
    258322!
    259323!-- Function result
    260     COMPLEX(wp), &
    261          DIMENSION(SIZE(array, 1), SIZE(array, 2), SIZE(array, 3)):: ft
    262 
    263     INTEGER(iwp) ::  ishape(3)
    264 
    265 !
    266 !-- Intrinsics used
    267     INTRINSIC SIZE, SHAPE
    268 
    269     ft = array
    270     ishape = SHAPE( array)
    271     CALL fftn(ft, ishape, dim, inv, stat)
    272 
    273  END FUNCTION fft3d
    274 
    275 
    276 !------------------------------------------------------------------------------!
    277 ! Description:
    278 ! ------------
    279 !> @todo Missing function description.
    280 !------------------------------------------------------------------------------!
    281  FUNCTION fft4d(array, dim, inv, stat) RESULT(ft)
    282 !
    283 !-- Formal parameters
    284     COMPLEX(wp),  DIMENSION(:,:,:,:), INTENT(IN)           :: array
    285     INTEGER(iwp), DIMENSION(:),       INTENT(IN),  OPTIONAL:: dim
    286     INTEGER(iwp),                     INTENT(OUT), OPTIONAL:: stat
    287     LOGICAL,                          INTENT(IN),  OPTIONAL:: inv
    288 !
    289 !-- Function result
    290     COMPLEX(wp), DIMENSION( &
    291          SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4)):: ft
    292 
    293     INTEGER(iwp) ::  ishape(4)
     324    COMPLEX(wp), DIMENSION(SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), SIZE(array, 5)) ::  ft  !<
     325
     326    INTEGER(iwp) ::  ishape(5)  !<
     327
    294328!
    295329!-- Intrinsics used
    … …  
    300334    CALL fftn(ft, ishape, dim, inv, stat)
    301335
    302  END FUNCTION fft4d
    303 
    304 
    305 !------------------------------------------------------------------------------!
     336 END FUNCTION fft5d
     337
     338
     339!--------------------------------------------------------------------------------------------------!
    306340! Description:
    307341! ------------
    308342!> @todo Missing function description.
    309 !------------------------------------------------------------------------------!
    310  FUNCTION fft5d(array, dim, inv, stat) RESULT(ft)
     343!--------------------------------------------------------------------------------------------------!
     344 FUNCTION fft6d( array, dim, inv, stat ) RESULT( ft )
    311345!
    312346!-- Formal parameters
    313     COMPLEX(wp),  DIMENSION(:,:,:,:,:), INTENT(IN)           :: array
    314     INTEGER(iwp), DIMENSION(:),         INTENT(IN),  OPTIONAL:: dim
    315     INTEGER(iwp),                       INTENT(OUT), OPTIONAL:: stat
    316     LOGICAL,                            INTENT(IN),  OPTIONAL:: inv
     347    COMPLEX(wp),  DIMENSION(:,:,:,:,:,:), INTENT(IN)            ::  array  !<
     348    INTEGER(iwp), DIMENSION(:),           INTENT(IN),  OPTIONAL ::  dim    !<
     349    INTEGER(iwp),                         INTENT(OUT), OPTIONAL ::  stat   !<
     350    LOGICAL,                              INTENT(IN),  OPTIONAL ::  inv    !<
    317351!
    318352!-- Function result
    319     COMPLEX(wp), DIMENSION( &
    320          SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), &
    321          SIZE(array, 5)):: ft
    322 
    323     INTEGER(iwp) ::  ishape(5)
     353    COMPLEX(wp), DIMENSION( SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4),        &
     354                            SIZE(array, 5), SIZE(array, 6)) ::  ft  !<
     355
     356    INTEGER(iwp) ::  ishape(6)  !<
    324357
    325358!
    … …  
    331364    CALL fftn(ft, ishape, dim, inv, stat)
    332365
    333  END FUNCTION fft5d
    334 
    335 
    336 !------------------------------------------------------------------------------!
     366 END FUNCTION fft6d
     367
     368
     369!--------------------------------------------------------------------------------------------------!
    337370! Description:
    338371! ------------
    339372!> @todo Missing function description.
    340 !------------------------------------------------------------------------------!
    341  FUNCTION fft6d(array, dim, inv, stat) RESULT(ft)
     373!--------------------------------------------------------------------------------------------------!
     374 FUNCTION fft7d( array, dim, inv, stat ) RESULT( ft )
    342375!
    343376!-- Formal parameters
    344     COMPLEX(wp),  DIMENSION(:,:,:,:,:,:), INTENT(IN)           :: array
    345     INTEGER(iwp), DIMENSION(:),           INTENT(IN),  OPTIONAL:: dim
    346     INTEGER(iwp),                         INTENT(OUT), OPTIONAL:: stat
    347     LOGICAL,                              INTENT(IN),  OPTIONAL:: inv
     377    COMPLEX(wp), DIMENSION(:,:,:,:,:,:,:), INTENT(IN)            ::  array  !<
     378    INTEGER(iwp),            DIMENSION(:), INTENT(IN),  OPTIONAL ::  dim    !<
     379    INTEGER(iwp),                          INTENT(OUT), OPTIONAL ::  stat   !<
     380    LOGICAL,                               INTENT(IN),  OPTIONAL ::  inv    !<
    348381!
    349382!-- Function result
    350     COMPLEX(wp), DIMENSION( &
    351          SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), &
    352          SIZE(array, 5), SIZE(array, 6)):: ft
    353 
    354     INTEGER(iwp) ::  ishape(6)
     383    COMPLEX(wp), DIMENSION( SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4),        &
     384                            SIZE(array, 5), SIZE(array, 6), SIZE(array, 7)) ::  ft  !<
     385
     386    INTEGER(iwp) ::  ishape(7)  !<
    355387
    356388!
    … …  
    362394    CALL fftn(ft, ishape, dim, inv, stat)
    363395
    364  END FUNCTION fft6d
    365 
    366 
    367 !------------------------------------------------------------------------------!
    368 ! Description:
    369 ! ------------
    370 !> @todo Missing function description.
    371 !------------------------------------------------------------------------------!
    372  FUNCTION fft7d(array, dim, inv, stat) RESULT(ft)
     396 END FUNCTION fft7d
     397
     398
     399!--------------------------------------------------------------------------------------------------!
     400! Description:
     401! ------------
     402!> @todo Missing subroutine description.
     403!--------------------------------------------------------------------------------------------------!
     404 SUBROUTINE fftn( array, shape, dim, inv, stat )
    373405!
    374406!-- Formal parameters
    375     COMPLEX(wp), DIMENSION(:,:,:,:,:,:,:), INTENT(IN)           :: array
    376     INTEGER(iwp),          DIMENSION(:),   INTENT(IN),  OPTIONAL:: dim
    377     INTEGER(iwp),                          INTENT(OUT), OPTIONAL:: stat
    378     LOGICAL,                               INTENT(IN),  OPTIONAL:: inv
    379 !
    380 !-- Function result
    381     COMPLEX(wp), DIMENSION( &
    382          SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), &
    383          SIZE(array, 5), SIZE(array, 6), SIZE(array, 7)):: ft
    384 
    385     INTEGER(iwp) ::  ishape(7)
    386 
    387 !
    388 !-- Intrinsics used
    389     INTRINSIC SIZE, SHAPE
    390 
    391     ft = array
    392     ishape = SHAPE( array )
    393     CALL fftn(ft, ishape, dim, inv, stat)
    394 
    395  END FUNCTION fft7d
    396 
    397 
    398 !------------------------------------------------------------------------------!
    399 ! Description:
    400 ! ------------
    401 !> @todo Missing subroutine description.
    402 !------------------------------------------------------------------------------!
    403  SUBROUTINE fftn(array, shape, dim, inv, stat)
    404 !
    405 !-- Formal parameters
    406     COMPLEX(wp),  DIMENSION(*), INTENT(INOUT)        :: array
    407     INTEGER(iwp), DIMENSION(:), INTENT(IN)           :: shape
    408     INTEGER(iwp), DIMENSION(:), INTENT(IN),  OPTIONAL:: dim
    409     INTEGER(iwp),               INTENT(OUT), OPTIONAL:: stat
    410     LOGICAL,                    INTENT(IN),  OPTIONAL:: inv
     407    COMPLEX(wp),  DIMENSION(*), INTENT(INOUT)         ::  array  !<
     408    INTEGER(iwp), DIMENSION(:), INTENT(IN)            ::  shape  !<
     409    INTEGER(iwp), DIMENSION(:), INTENT(IN),  OPTIONAL ::  dim    !<
     410    INTEGER(iwp),               INTENT(OUT), OPTIONAL ::  stat   !<
     411    LOGICAL,                    INTENT(IN),  OPTIONAL ::  inv    !<
    411412!
    412413!-- Local arrays
    413     INTEGER(iwp), DIMENSION(SIZE(shape)):: d
     414    INTEGER(iwp), DIMENSION(SIZE(shape)) ::  d  !<
    414415!
    415416!-- Local scalars
    416     LOGICAL      :: inverse
    417     INTEGER(iwp) :: i, ndim, ntotal
    418     REAL(wp):: scale
     417    LOGICAL      ::  inverse          !<
     418    INTEGER(iwp) ::  i, ndim, ntotal  !<
     419    REAL(wp)     ::  scale            !<
    419420!
    420421!-- Intrinsics used
    … …  
    423424!
    424425!-- Optional parameter settings
    425     IF (PRESENT(inv)) THEN
     426    IF ( PRESENT( inv ) ) THEN
    426427       inverse = inv
    427428    ELSE
    428429       inverse = .FALSE.
    429430    END IF
    430     IF (PRESENT(dim)) THEN
    431        ndim = MIN(SIZE(dim), SIZE(d))
    432        d(1:ndim) = DIM(1:ndim)
     431    IF ( PRESENT( dim ) ) THEN
     432       ndim = MIN( SIZE( dim ), SIZE( d ) )
     433       d(1:ndim) = DIM( 1:ndim )
    433434    ELSE
    434        ndim = SIZE(d)
    435        d = (/(i, i = 1, SIZE(d))/)
     435       ndim = SIZE( d )
     436        d = (/( i, i = 1, SIZE( d ) )/)
    436437    END IF
    437438
    438     ntotal = PRODUCT(shape)
    439     scale = SQRT(1.0_wp / PRODUCT(shape(d(1:ndim))))
    440     DO i = 1, ntotal
    441        array(i) = CMPLX(REAL(array(i)) * scale, AIMAG(array(i)) * scale, &
    442             KIND=wp)
     439    ntotal = PRODUCT( shape )
     440    scale = SQRT( 1.0_wp / PRODUCT( shape( d(1:ndim) ) ) )
     441    DO  i = 1, ntotal
     442       array(i) = CMPLX( REAL( array(i) ) * scale, AIMAG( array(i) ) * scale, KIND = wp )
    443443    END DO
    444444
    445     DO i = 1, ndim
    446        CALL fftradix(array, ntotal, shape(d(i)), PRODUCT(shape(1:d(i))), &
    447             inverse, stat)
    448        IF (PRESENT(stat)) THEN
    449           IF (stat /=0) RETURN
     445    DO  i = 1, ndim
     446       CALL fftradix( array, ntotal, shape( d(i) ), PRODUCT( shape( 1:d(i) ) ), inverse, stat )
     447       IF ( PRESENT( stat ) )  THEN
     448          IF ( stat /= 0 )  RETURN
    450449       END IF
    451450    END DO
    … …  
    454453
    455454
    456 !------------------------------------------------------------------------------!
     455!--------------------------------------------------------------------------------------------------!
    457456! Description:
    458457! ------------
    459458!> @todo Missing subroutine description.
    460 !------------------------------------------------------------------------------!
    461  SUBROUTINE fftradix(array, ntotal, npass, nspan, inv, stat)
     459!--------------------------------------------------------------------------------------------------!
     460 SUBROUTINE fftradix( array, ntotal, npass, nspan, inv, stat )
    462461!
    463462!-- Formal parameters
    464     COMPLEX(wp),  DIMENSION(*), INTENT(INOUT)        :: array
    465     INTEGER(iwp),               INTENT(IN)           :: ntotal, npass, nspan
    466     INTEGER(iwp),               INTENT(OUT), OPTIONAL:: stat
    467     LOGICAL,                    INTENT(IN)           :: inv
     463    COMPLEX(wp),  DIMENSION(*), INTENT(INOUT)         ::  array                 !<
     464    INTEGER(iwp),               INTENT(IN)            ::  ntotal, npass, nspan  !<
     465    INTEGER(iwp),               INTENT(OUT), OPTIONAL ::  stat                  !<
     466    LOGICAL,                    INTENT(IN)            ::  inv                   !<
    468467!
    469468!-- Local arrays
    470     COMPLEX(wp),  DIMENSION(:), ALLOCATABLE  :: ctmp
    471     INTEGER(iwp), DIMENSION(BIT_SIZE(0))     :: factor
    472     INTEGER(iwp), DIMENSION(:), ALLOCATABLE  :: perm
    473     REAL(wp),     DIMENSION(:), ALLOCATABLE  :: sine, cosine
     469    COMPLEX(wp),  DIMENSION(:), ALLOCATABLE ::  ctmp          !<
     470    INTEGER(iwp), DIMENSION(BIT_SIZE(0))    ::  factor        !<
     471    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  perm          !<
     472    REAL(wp),     DIMENSION(:), ALLOCATABLE ::  sine, cosine  !<
    474473!
    475474!-- Local scalars
    476     INTEGER(iwp)         :: maxfactor, nfactor, nsquare, nperm
     475    INTEGER(iwp) ::  maxfactor, nfactor, nsquare, nperm  !<
    477476!
    478477!-- Intrinsics used
    479     INTRINSIC MAXVAL, MOD, PRESENT, ISHFT, BIT_SIZE, SIN, COS, &
    480          CMPLX, REAL, AIMAG
    481 
    482     IF (npass <= 1) RETURN
    483 
    484     CALL factorize(npass, factor, nfactor, nsquare)
    485 
    486     maxfactor = MAXVAL(factor(:nfactor))
    487     IF (nfactor - ISHFT(nsquare, 1) > 0) THEN
    488        nperm = MAX(nfactor + 1, PRODUCT(factor(nsquare+1: nfactor-nsquare)) - 1)
     478    INTRINSIC MAXVAL, MOD, PRESENT, ISHFT, BIT_SIZE, SIN, COS, CMPLX, REAL, AIMAG
     479
     480    IF ( npass <= 1 )  RETURN
     481
     482    CALL factorize( npass, factor, nfactor, nsquare )
     483
     484    maxfactor = MAXVAL( factor(:nfactor) )
     485    IF ( nfactor - ISHFT( nsquare, 1 ) > 0 )  THEN
     486       nperm = MAX( nfactor + 1, PRODUCT( factor(nsquare+1: nfactor-nsquare) ) - 1 )
    489487    ELSE
    490488       nperm = nfactor + 1
    491489    END IF
    492490
    493     IF (PRESENT(stat)) THEN
    494        ALLOCATE(ctmp(maxfactor), sine(maxfactor), cosine(maxfactor), STAT=stat)
    495        IF (stat /= 0) RETURN
    496        CALL transform(array, ntotal, npass, nspan, &
    497             factor, nfactor, ctmp, sine, cosine, inv)
    498        DEALLOCATE(sine, cosine, STAT=stat)
    499        IF (stat /= 0) RETURN
    500        ALLOCATE(perm(nperm), STAT=stat)
    501        IF (stat /= 0) RETURN
    502        CALL permute(array, ntotal, npass, nspan, &
    503             factor, nfactor, nsquare, maxfactor, &
    504             ctmp, perm)
    505        DEALLOCATE(perm, ctmp, STAT=stat)
    506        IF (stat /= 0) RETURN
     491    IF ( PRESENT( stat ) )  THEN
     492       ALLOCATE( ctmp(maxfactor), sine(maxfactor), cosine(maxfactor), STAT = stat )
     493       IF ( stat /= 0 )  RETURN
     494       CALL transform( array, ntotal, npass, nspan, factor, nfactor, ctmp, sine, cosine, inv )
     495       DEALLOCATE( sine, cosine, STAT = stat )
     496       IF ( stat /= 0 )  RETURN
     497       ALLOCATE( perm(nperm), STAT = stat )
     498       IF ( stat /= 0 )  RETURN
     499       CALL permute( array, ntotal, npass, nspan, factor, nfactor, nsquare, maxfactor, ctmp, perm )
     500       DEALLOCATE( perm, ctmp, STAT = stat )
     501       IF ( stat /= 0 )  RETURN
    507502    ELSE
    508        ALLOCATE(ctmp(maxfactor), sine(maxfactor), cosine(maxfactor))
    509        CALL transform(array, ntotal, npass, nspan, &
    510             factor, nfactor, ctmp, sine, cosine, inv)
    511        DEALLOCATE(sine, cosine)
    512        ALLOCATE(perm(nperm))
    513        CALL permute(array, ntotal, npass, nspan, &
    514             factor, nfactor, nsquare, maxfactor, &
    515             ctmp, perm)
    516        DEALLOCATE(perm, ctmp)
     503       ALLOCATE( ctmp(maxfactor), sine(maxfactor), cosine(maxfactor) )
     504       CALL transform( array, ntotal, npass, nspan, factor, nfactor, ctmp, sine, cosine, inv )
     505       DEALLOCATE( sine, cosine )
     506       ALLOCATE( perm(nperm) )
     507       CALL permute( array, ntotal, npass, nspan, factor, nfactor, nsquare, maxfactor, ctmp, perm )
     508       DEALLOCATE( perm, ctmp )
    517509    END IF
    518510
    519511
    520   CONTAINS
    521 
    522 
    523 !------------------------------------------------------------------------------!
     512 CONTAINS
     513
     514
     515!--------------------------------------------------------------------------------------------------!
    524516! Description:
    525517! ------------
    526518!> @todo Missing subroutine description.
    527 !------------------------------------------------------------------------------!
    528     SUBROUTINE factorize(npass, factor, nfactor, nsquare)
    529 !
    530 !--   Formal parameters
    531       INTEGER(iwp),               INTENT(IN) :: npass
    532       INTEGER(iwp), DIMENSION(*), INTENT(OUT):: factor
    533       INTEGER(iwp),               INTENT(OUT):: nfactor, nsquare
    534 !
    535 !--   Local scalars
    536       INTEGER(iwp):: j, jj, k
    537 
    538       nfactor = 0
    539       k = npass
    540       DO WHILE (MOD(k, 16) == 0)
    541          nfactor = nfactor + 1
    542          factor(nfactor) = 4
    543          k = k / 16
    544       END DO
    545       j = 3
    546       jj = 9
    547       DO
    548          DO WHILE (MOD(k, jj) == 0)
    549             nfactor = nfactor + 1
    550             factor(nfactor) = j
    551             k = k / jj
    552          END DO
    553          j = j + 2
    554          jj = j * j
    555          IF (jj > k) EXIT
    556       END DO
    557       IF (k <= 4) THEN
    558          nsquare = nfactor
    559          factor(nfactor + 1) = k
    560          IF (k /= 1) nfactor = nfactor + 1
    561       ELSE
    562          IF (k - ISHFT(k / 4, 2) == 0) THEN
    563             nfactor = nfactor + 1
    564             factor(nfactor) = 2
    565             k = k / 4
    566          END IF
    567          nsquare = nfactor
    568          j = 2
    569          DO
    570             IF (MOD(k, j) == 0) THEN
    571                nfactor = nfactor + 1
    572                factor(nfactor) = j
    573                k = k / j
    574             END IF
    575             j = ISHFT((j + 1) / 2, 1) + 1
    576             IF (j > k) EXIT
    577          END DO
    578       END IF
    579       IF (nsquare > 0) THEN
    580          j = nsquare
    581          DO
    582             nfactor = nfactor + 1
    583             factor(nfactor) = factor(j)
    584             j = j - 1
    585             IF (j==0) EXIT
    586          END DO
    587       END IF
    588 
    589     END SUBROUTINE factorize
    590 
    591 
    592 !------------------------------------------------------------------------------!
     519!--------------------------------------------------------------------------------------------------!
     520 SUBROUTINE factorize( npass, factor, nfactor, nsquare )
     521!
     522!-- Formal parameters
     523    INTEGER(iwp),               INTENT(IN)  ::  npass             !<
     524    INTEGER(iwp),               INTENT(OUT) ::  nfactor, nsquare  !<
     525    INTEGER(iwp), DIMENSION(*), INTENT(OUT) ::  factor            !<
     526!
     527!-- Local scalars
     528    INTEGER(iwp) ::  j, jj, k  !<
     529
     530    nfactor = 0
     531    k = npass
     532    DO WHILE ( MOD( k, 16 ) == 0 )
     533       nfactor = nfactor + 1
     534       factor(nfactor) = 4
     535       k = k / 16
     536    END DO
     537    j = 3
     538    jj = 9
     539    DO
     540       DO WHILE ( MOD( k, jj ) == 0 )
     541          nfactor = nfactor + 1
     542          factor(nfactor) = j
     543          k = k / jj
     544       END DO
     545       j = j + 2
     546       jj = j * j
     547       IF ( jj > k ) EXIT
     548    END DO
     549    IF ( k <= 4 ) THEN
     550       nsquare = nfactor
     551       factor(nfactor + 1) = k
     552       IF ( k /= 1 ) nfactor = nfactor + 1
     553    ELSE
     554       IF ( k - ISHFT( k / 4, 2 ) == 0 ) THEN
     555          nfactor = nfactor + 1
     556          factor(nfactor) = 2
     557          k = k / 4
     558       END IF
     559       nsquare = nfactor
     560       j = 2
     561       DO
     562          IF ( MOD(k, j) == 0 ) THEN
     563             nfactor = nfactor + 1
     564             factor(nfactor) = j
     565             k = k / j
     566          END IF
     567          j = ISHFT( (j + 1) / 2, 1 ) + 1
     568          IF ( j > k ) EXIT
     569       END DO
     570    END IF
     571    IF ( nsquare > 0 ) THEN
     572       j = nsquare
     573       DO
     574          nfactor = nfactor + 1
     575          factor(nfactor) = factor(j)
     576          j = j - 1
     577          IF ( j == 0 ) EXIT
     578       END DO
     579    END IF
     580
     581 END SUBROUTINE factorize
     582
     583
     584!--------------------------------------------------------------------------------------------------!
    593585! Description:
    594586! ------------
    595587!> @todo Missing subroutine description.
    596 !------------------------------------------------------------------------------!
    597     SUBROUTINE transform(array, ntotal, npass, nspan, &
    598          factor, nfactor, ctmp, sine, cosine, inv) !-- compute fourier transform
    599 !
    600 !--   Formal parameters
    601       COMPLEX(wp),  DIMENSION(*), INTENT(IN OUT):: array
    602       COMPLEX(wp),  DIMENSION(*), INTENT(OUT)   :: ctmp
    603       INTEGER(iwp),               INTENT(IN)    :: ntotal, npass, nspan
    604       INTEGER(iwp), DIMENSION(*), INTENT(IN)    :: factor
    605       INTEGER(iwp),               INTENT(IN)    :: nfactor
    606       LOGICAL,                    INTENT(IN)    :: inv
    607       REAL(wp),     DIMENSION(*), INTENT(OUT)   :: sine, cosine
    608 !
    609 !--   Local scalars
    610       INTEGER(iwp):: ii, ispan
    611       INTEGER(iwp):: j, jc, jf, jj
    612       INTEGER(iwp):: k, kk, kspan, k1, k2, k3, k4
    613       INTEGER(iwp):: nn, nt
    614       REAL(wp)    :: s60, c72, s72, pi2, radf
    615       REAL(wp)    :: c1, s1, c2, s2, c3, s3, cd, sd, ak
    616       COMPLEX(wp) :: cc, cj, ck, cjp, cjm, ckp, ckm
    617 
    618       c72 = cos72
    619       IF (inv) THEN
    620          s72 = sin72
    621          s60 = sin60
    622          pi2 = pi
    623       ELSE
    624          s72 = -sin72
    625          s60 = -sin60
    626          pi2 = -pi
    627       END IF
    628 
    629       nt = ntotal
    630       nn = nt - 1
    631       kspan = nspan
    632       jc = nspan / npass
    633       radf = pi2 * jc
    634       pi2 = pi2 * 2.0_wp !-- use 2 PI from here on
    635 
    636       ii = 0
    637       jf = 0
    638       DO
    639          sd = radf / kspan
    640          cd = SIN(sd)
    641          cd = 2.0_wp * cd * cd
    642          sd = SIN(sd + sd)
    643          kk = 1
    644          ii = ii + 1
    645 
    646          SELECT CASE (factor(ii))
    647          CASE (2)
    648 !
    649 !--         Transform for factor of 2 (including rotation factor)
    650             kspan = kspan / 2
    651             k1 = kspan + 2
    652             DO
    653                DO
    654                   k2 = kk + kspan
    655                   ck = array(k2)
    656                   array(k2) = array(kk)-ck
    657                   array(kk) = array(kk) + ck
    658                   kk = k2 + kspan
    659                   IF (kk > nn) EXIT
    660                END DO
    661                kk = kk - nn
    662                IF (kk > jc) EXIT
    663             END DO
    664             IF (kk > kspan) RETURN
    665             DO
    666                c1 = 1.0_wp - cd
    667                s1 = sd
    668                DO
    669                   DO
    670                      DO
    671                         k2 = kk + kspan
    672                         ck = array(kk) - array(k2)
    673                         array(kk) = array(kk) + array(k2)
    674                         array(k2) = ck * CMPLX(c1, s1, KIND=wp)
    675                         kk = k2 + kspan
    676                         IF (kk >= nt) EXIT
    677                      END DO
    678                      k2 = kk - nt
    679                      c1 = -c1
    680                      kk = k1 - k2
    681                      IF (kk <= k2) EXIT
    682                   END DO
    683                   ak = c1 - (cd * c1 + sd * s1)
    684                   s1 = sd * c1 - cd * s1 + s1
    685                   c1 = 2.0_wp - (ak * ak + s1 * s1)
    686                   s1 = s1 * c1
    687                   c1 = c1 * ak
    688                   kk = kk + jc
    689                   IF (kk >= k2) EXIT
    690                END DO
    691                k1 = k1 + 1 + 1
    692                kk = (k1 - kspan) / 2 + jc
    693                IF (kk > jc + jc) EXIT
    694             END DO
    695 
    696          CASE (4) !-- transform for factor of 4
    697             ispan = kspan
    698             kspan = kspan / 4
    699 
    700             DO
    701                c1 = 1.0_wp
    702                s1 = 0.0_wp
    703                DO
    704                   DO
    705                      k1 = kk + kspan
    706                      k2 = k1 + kspan
    707                      k3 = k2 + kspan
    708                      ckp = array(kk) + array(k2)
    709                      ckm = array(kk) - array(k2)
    710                      cjp = array(k1) + array(k3)
    711                      cjm = array(k1) - array(k3)
    712                      array(kk) = ckp + cjp
    713                      cjp = ckp - cjp
    714                      IF (inv) THEN
    715                         ckp = ckm + CMPLX(-AIMAG(cjm), REAL(cjm), KIND=wp)
    716                         ckm = ckm + CMPLX(AIMAG(cjm), -REAL(cjm), KIND=wp)
    717                      ELSE
    718                         ckp = ckm + CMPLX(AIMAG(cjm), -REAL(cjm), KIND=wp)
    719                         ckm = ckm + CMPLX(-AIMAG(cjm), REAL(cjm), KIND=wp)
    720                      END IF
    721 !
    722 !--                  Avoid useless multiplies
    723                      IF (s1 == 0.0_wp) THEN
    724                         array(k1) = ckp
    725                         array(k2) = cjp
    726                         array(k3) = ckm
    727                      ELSE
    728                         array(k1) = ckp * CMPLX(c1, s1, KIND=wp)
    729                         array(k2) = cjp * CMPLX(c2, s2, KIND=wp)
    730                         array(k3) = ckm * CMPLX(c3, s3, KIND=wp)
    731                      END IF
    732                      kk = k3 + kspan
    733                      IF (kk > nt) EXIT
    734                   END DO
    735 
    736                   c2 = c1 - (cd * c1 + sd * s1)
    737                   s1 = sd * c1 - cd * s1 + s1
    738                   c1 = 2.0_wp - (c2 * c2 + s1 * s1)
    739                   s1 = s1 * c1
    740                   c1 = c1 * c2
    741 !
    742 !--               Values of c2, c3, s2, s3 that will get used next time
    743                   c2 = c1 * c1 - s1 * s1
    744                   s2 = 2.0_wp * c1 * s1
    745                   c3 = c2 * c1 - s2 * s1
    746                   s3 = c2 * s1 + s2 * c1
    747                   kk = kk - nt + jc
    748                   IF (kk > kspan) EXIT
    749                END DO
    750                kk = kk - kspan + 1
    751                IF (kk > jc) EXIT
    752             END DO
    753             IF (kspan == jc) RETURN
    754 
    755          CASE default
    756 !
    757 !--         Transform for odd factors
    758             k = factor(ii)
    759             ispan = kspan
    760             kspan = kspan / k
    761 
    762             SELECT CASE (k)
    763             CASE (3) !-- transform for factor of 3 (optional code)
    764                DO
    765                   DO
    766                      k1 = kk + kspan
    767                      k2 = k1 + kspan
    768                      ck = array(kk)
    769                      cj = array(k1) + array(k2)
    770                      array(kk) = ck + cj
    771                      ck = ck - CMPLX( &
    772                           0.5_wp * REAL (cj), &
    773                           0.5_wp * AIMAG(cj), &
    774                           KIND=wp)
    775                      cj = CMPLX( &
    776                           (REAL (array(k1)) - REAL (array(k2))) * s60, &
    777                           (AIMAG(array(k1)) - AIMAG(array(k2))) * s60, &
    778                           KIND=wp)
    779                      array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND=wp)
    780                      array(k2) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND=wp)
    781                      kk = k2 + kspan
    782                      IF (kk >= nn) EXIT
    783                   END DO
    784                   kk = kk - nn
    785                   IF (kk > kspan) EXIT
    786                END DO
    787 
    788             CASE (5) !-- transform for factor of 5 (optional code)
    789                c2 = c72 * c72 - s72 * s72
    790                s2 = 2.0_wp * c72 * s72
    791                DO
    792                   DO
    793                      k1 = kk + kspan
    794                      k2 = k1 + kspan
    795                      k3 = k2 + kspan
    796                      k4 = k3 + kspan
    797                      ckp = array(k1) + array(k4)
    798                      ckm = array(k1) - array(k4)
    799                      cjp = array(k2) + array(k3)
    800                      cjm = array(k2) - array(k3)
    801                      cc = array(kk)
    802                      array(kk) = cc + ckp + cjp
    803                      ck = CMPLX(REAL(ckp) * c72, AIMAG(ckp) * c72, &
    804                           KIND=wp) + &
    805                           CMPLX(REAL(cjp) * c2,  AIMAG(cjp) * c2,  &
    806                           KIND=wp) + cc
    807                      cj = CMPLX(REAL(ckm) * s72, AIMAG(ckm) * s72, &
    808                           KIND=wp) + &
    809                           CMPLX(REAL(cjm) * s2,  AIMAG(cjm) * s2,  &
    810                           KIND=wp)
    811                      array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND=wp)
    812                      array(k4) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND=wp)
    813                      ck = CMPLX(REAL(ckp) * c2,  AIMAG(ckp) * c2,  &
    814                           KIND=wp) + &
    815                           CMPLX(REAL(cjp) * c72, AIMAG(cjp) * c72, &
    816                           KIND=wp) + cc
    817                      cj = CMPLX(REAL(ckm) * s2,  AIMAG(ckm) * s2,  &
    818                           KIND=wp) - &
    819                           CMPLX(REAL(cjm) * s72, AIMAG(cjm) * s72, &
    820                           KIND=wp)
    821                      array(k2) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND=wp)
    822                      array(k3) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND=wp)
    823                      kk = k4 + kspan
    824                      IF (kk >= nn) EXIT
    825                   END DO
    826                   kk = kk - nn
    827                   IF (kk > kspan) EXIT
    828                END DO
    829 
    830             CASE default
    831                IF (k /= jf) THEN
    832                   jf = k
    833                   s1 = pi2 / k
    834                   c1 = COS(s1)
    835                   s1 = SIN(s1)
    836                   cosine (jf) = 1.0_wp
    837                   sine (jf) = 0.0_wp
    838                   j = 1
    839                   DO
    840                      cosine (j) = cosine (k) * c1 + sine (k) * s1
    841                      sine (j) = cosine (k) * s1 - sine (k) * c1
    842                      k = k-1
    843                      cosine (k) = cosine (j)
    844                      sine (k) = -sine (j)
    845                      j = j + 1
    846                      IF (j >= k) EXIT
    847                   END DO
    848                END IF
    849                DO
    850                   DO
    851                      k1 = kk
    852                      k2 = kk + ispan
    853                      cc = array(kk)
    854                      ck = cc
    855                      j = 1
    856                      k1 = k1 + kspan
    857                      DO
    858                         k2 = k2 - kspan
    859                         j = j + 1
    860                         ctmp(j) = array(k1) + array(k2)
    861                         ck = ck + ctmp(j)
    862                         j = j + 1
    863                         ctmp(j) = array(k1) - array(k2)
    864                         k1 = k1 + kspan
    865                         IF (k1 >= k2) EXIT
    866                      END DO
    867                      array(kk) = ck
    868                      k1 = kk
    869                      k2 = kk + ispan
    870                      j = 1
    871                      DO
    872                         k1 = k1 + kspan
    873                         k2 = k2 - kspan
    874                         jj = j
    875                         ck = cc
    876                         cj = (0.0_wp, 0.0_wp)
    877                         k = 1
    878                         DO
    879                            k = k + 1
    880                            ck = ck + CMPLX( &
    881                                 REAL (ctmp(k)) * cosine(jj), &
    882                                 AIMAG(ctmp(k)) * cosine(jj), KIND=wp)
    883                            k = k + 1
    884                            cj = cj + CMPLX( &
    885                                 REAL (ctmp(k)) * sine(jj), &
    886                                 AIMAG(ctmp(k)) * sine(jj), KIND=wp)
    887                            jj = jj + j
    888                            IF (jj > jf) jj = jj - jf
    889                            IF (k >= jf) EXIT
    890                         END DO
    891                         k = jf - j
    892                         array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), &
    893                              KIND=wp)
    894                         array(k2) = ck + CMPLX(AIMAG(cj), -REAL(cj), &
    895                              KIND=wp)
    896                         j = j + 1
    897                         IF (j >= k) EXIT
    898                      END DO
    899                      kk = kk + ispan
    900                      IF (kk > nn) EXIT
    901                   END DO
    902                   kk = kk - nn
    903                   IF (kk > kspan) EXIT
    904                END DO
    905 
    906             END SELECT
    907 !
    908 !--         Multiply by rotation factor (except for factors of 2 and 4)
    909             IF (ii == nfactor) RETURN
    910             kk = jc + 1
    911             DO
    912                c2 = 1.0_wp - cd
    913                s1 = sd
    914                DO
    915                   c1 = c2
    916                   s2 = s1
    917                   kk = kk + kspan
    918                   DO
    919                      DO
    920                         array(kk) = CMPLX(c2, s2, KIND=wp) * array(kk)
    921                         kk = kk + ispan
    922                         IF (kk > nt) EXIT
    923                      END DO
    924                      ak = s1 * s2
    925                      s2 = s1 * c2 + c1 * s2
    926                      c2 = c1 * c2 - ak
    927                      kk = kk - nt + kspan
    928                      IF (kk > ispan) EXIT
    929                   END DO
    930                   c2 = c1 - (cd * c1 + sd * s1)
    931                   s1 = s1 + sd * c1 - cd * s1
    932                   c1 = 2.0_wp - (c2 * c2 + s1 * s1)
    933                   s1 = s1 * c1
    934                   c2 = c2 * c1
    935                   kk = kk - ispan + jc
    936                   IF (kk > kspan) EXIT
    937                END DO
    938                kk = kk - kspan + jc + 1
    939                IF (kk > jc + jc) EXIT
    940             END DO
    941 
    942          END SELECT
    943       END DO
    944     END SUBROUTINE transform
    945 
    946 
    947 !------------------------------------------------------------------------------!
     588!--------------------------------------------------------------------------------------------------!
     589 SUBROUTINE transform( array, ntotal, npass, nspan, factor, nfactor, ctmp, sine, cosine, inv )
     590!-- Compute fourier transform
     591
     592!
     593!-- Formal parameters
     594    COMPLEX(wp),  DIMENSION(*), INTENT(IN OUT) ::  array                 !<
     595    COMPLEX(wp),  DIMENSION(*), INTENT(OUT)    ::  ctmp                  !<
     596    INTEGER(iwp),               INTENT(IN)     ::  ntotal, npass, nspan  !<
     597    INTEGER(iwp), DIMENSION(*), INTENT(IN)     ::  factor                !<
     598    INTEGER(iwp),               INTENT(IN)     ::  nfactor               !<
     599    LOGICAL,                    INTENT(IN)     ::  inv                   !<
     600    REAL(wp),     DIMENSION(*), INTENT(OUT)    ::  sine, cosine          !<
     601!
     602!-- Local scalars
     603    COMPLEX(wp)  ::  cc, cj, ck, cjp, cjm, ckp, ckm      !<
     604    INTEGER(iwp) ::  ii, ispan                           !<
     605    INTEGER(iwp) ::  j, jc, jf, jj                       !<
     606    INTEGER(iwp) ::  k, kk, kspan, k1, k2, k3, k4        !<
     607    INTEGER(iwp) ::  nn, nt                              !<
     608    REAL(wp)     ::  s60, c72, s72, pi2, radf            !<
     609    REAL(wp)     ::  c1, s1, c2, s2, c3, s3, cd, sd, ak  !<
     610
     611    c72 = cos72
     612    IF ( inv )  THEN
     613       s72 = sin72
     614       s60 = sin60
     615       pi2 = pi
     616    ELSE
     617       s72 = - sin72
     618       s60 = - sin60
     619       pi2 = - pi
     620    END IF
     621
     622    nt = ntotal
     623    nn = nt - 1
     624    kspan = nspan
     625    jc = nspan / npass
     626    radf = pi2 * jc
     627    pi2 = pi2 * 2.0_wp !-- Use 2 PI from here on
     628
     629    ii = 0
     630    jf = 0
     631    DO
     632       sd = radf / kspan
     633       cd = SIN( sd )
     634       cd = 2.0_wp * cd * cd
     635       sd = SIN( sd + sd )
     636       kk = 1
     637       ii = ii + 1
     638
     639       SELECT CASE ( factor(ii) )
     640       CASE ( 2 )
     641!
     642!--       Transform for factor of 2 (including rotation factor)
     643          kspan = kspan / 2
     644          k1 = kspan + 2
     645          DO
     646             DO
     647                k2 = kk + kspan
     648                ck = array(k2)
     649                array(k2) = array(kk) - ck
     650                array(kk) = array(kk) + ck
     651                kk = k2 + kspan
     652                IF ( kk > nn )  EXIT
     653             END DO
     654             kk = kk - nn
     655             IF ( kk > jc )  EXIT
     656          END DO
     657          IF ( kk > kspan )  RETURN
     658          DO
     659             c1 = 1.0_wp - cd
     660             s1 = sd
     661             DO
     662                DO
     663                   DO
     664                      k2 = kk + kspan
     665                      ck = array(kk) - array(k2)
     666                      array(kk) = array(kk) + array(k2)
     667                      array(k2) = ck * CMPLX( c1, s1, KIND = wp )
     668                      kk = k2 + kspan
     669                      IF ( kk >= nt )  EXIT
     670                   END DO
     671                   k2 = kk - nt
     672                   c1 = - c1
     673                   kk = k1 - k2
     674                   IF ( kk <= k2 )  EXIT
     675                END DO
     676                ak = c1 - (cd * c1 + sd * s1)
     677                s1 = sd * c1 - cd * s1 + s1
     678                c1 = 2.0_wp - ( ak * ak + s1 * s1 )
     679                s1 = s1 * c1
     680                c1 = c1 * ak
     681                kk = kk + jc
     682                IF ( kk >= k2 ) EXIT
     683             END DO
     684             k1 = k1 + 1 + 1
     685             kk = ( k1 - kspan ) / 2 + jc
     686             IF ( kk > jc + jc )  EXIT
     687          END DO
     688!
     689!--    Transform for factor of 4
     690       CASE ( 4 )
     691          ispan = kspan
     692          kspan = kspan / 4
     693
     694          DO
     695             c1 = 1.0_wp
     696             s1 = 0.0_wp
     697             DO
     698                DO
     699                   k1 = kk + kspan
     700                   k2 = k1 + kspan
     701                   k3 = k2 + kspan
     702                   ckp = array(kk) + array(k2)
     703                   ckm = array(kk) - array(k2)
     704                   cjp = array(k1) + array(k3)
     705                   cjm = array(k1) - array(k3)
     706                   array(kk) = ckp + cjp
     707                   cjp = ckp - cjp
     708                   IF ( inv )  THEN
     709                      ckp = ckm + CMPLX( - AIMAG( cjm ), REAL( cjm ), KIND = wp )
     710                      ckm = ckm + CMPLX( AIMAG( cjm ), - REAL( cjm ), KIND = wp )
     711                   ELSE
     712                      ckp = ckm + CMPLX( AIMAG( cjm ), - REAL( cjm ), KIND = wp )
     713                      ckm = ckm + CMPLX( - AIMAG( cjm ), REAL( cjm ), KIND = wp )
     714                   END IF
     715!
     716!--                Avoid useless multiplies
     717                   IF ( s1 == 0.0_wp )  THEN
     718                      array(k1) = ckp
     719                      array(k2) = cjp
     720                      array(k3) = ckm
     721                   ELSE
     722                      array(k1) = ckp * CMPLX( c1, s1, KIND = wp )
     723                      array(k2) = cjp * CMPLX( c2, s2, KIND = wp )
     724                      array(k3) = ckm * CMPLX( c3, s3, KIND = wp )
     725                   END IF
     726                   kk = k3 + kspan
     727                   IF ( kk > nt )  EXIT
     728                END DO
     729
     730                c2 = c1 - ( cd * c1 + sd * s1 )
     731                s1 = sd * c1 - cd * s1 + s1
     732                c1 = 2.0_wp - ( c2 * c2 + s1 * s1 )
     733                s1 = s1 * c1
     734                c1 = c1 * c2
     735!
     736!--             Values of c2, c3, s2, s3 that will get used next time
     737                c2 = c1 * c1 - s1 * s1
     738                s2 = 2.0_wp * c1 * s1
     739                c3 = c2 * c1 - s2 * s1
     740                s3 = c2 * s1 + s2 * c1
     741                kk = kk - nt + jc
     742                IF ( kk > kspan )  EXIT
     743             END DO
     744             kk = kk - kspan + 1
     745             IF ( kk > jc )  EXIT
     746          END DO
     747          IF ( kspan == jc )  RETURN
     748
     749       CASE default
     750!
     751!--       Transform for odd factors
     752          k = factor(ii)
     753          ispan = kspan
     754          kspan = kspan / k
     755
     756          SELECT CASE ( k )
     757!
     758!--       Transform for factor of 3 (optional code)
     759          CASE ( 3 )
     760             DO
     761                DO
     762                   k1 = kk + kspan
     763                   k2 = k1 + kspan
     764                   ck = array(kk)
     765                   cj = array(k1) + array(k2)
     766                   array(kk) = ck + cj
     767                   ck = ck - CMPLX( 0.5_wp * REAL( cj ), 0.5_wp * AIMAG( cj ), KIND = wp )
     768                   cj = CMPLX( ( REAL( array(k1) ) - REAL( array(k2) ) ) * s60,                    &
     769                               ( AIMAG( array(k1) ) - AIMAG( array(k2) ) ) * s60, KIND = wp )
     770                   array(k1) = ck + CMPLX( - AIMAG( cj ), REAL( cj ), KIND = wp )
     771                   array(k2) = ck + CMPLX( AIMAG( cj ), - REAL( cj ), KIND = wp )
     772                   kk = k2 + kspan
     773                   IF ( kk >= nn )  EXIT
     774                END DO
     775                kk = kk - nn
     776                IF ( kk > kspan )  EXIT
     777             END DO
     778!
     779!--       Transform for factor of 5 (optional code)
     780          CASE ( 5 )
     781             c2 = c72 * c72 - s72 * s72
     782             s2 = 2.0_wp * c72 * s72
     783             DO
     784                DO
     785                   k1 = kk + kspan
     786                   k2 = k1 + kspan
     787                   k3 = k2 + kspan
     788                   k4 = k3 + kspan
     789                   ckp = array(k1) + array(k4)
     790                   ckm = array(k1) - array(k4)
     791                   cjp = array(k2) + array(k3)
     792                   cjm = array(k2) - array(k3)
     793                   cc = array(kk)
     794                   array(kk) = cc + ckp + cjp
     795                   ck = CMPLX( REAL( ckp ) * c72, AIMAG( ckp ) * c72, KIND = wp ) +                &
     796                        CMPLX( REAL( cjp ) * c2,  AIMAG( cjp ) * c2,  KIND = wp ) + cc
     797                   cj = CMPLX( REAL( ckm ) * s72, AIMAG( ckm ) * s72, KIND = wp) +                 &
     798                        CMPLX( REAL( cjm ) * s2,  AIMAG( cjm ) * s2,  KIND = wp )
     799                   array(k1) = ck + CMPLX( - AIMAG( cj ), REAL( cj ), KIND = wp )
     800                   array(k4) = ck + CMPLX( AIMAG( cj ), - REAL( cj ), KIND = wp )
     801                   ck = CMPLX( REAL( ckp ) * c2,  AIMAG( ckp ) * c2,  KIND = wp ) +                &
     802                        CMPLX( REAL( cjp ) * c72, AIMAG( cjp ) * c72, KIND = wp ) + cc
     803                   cj = CMPLX( REAL( ckm ) * s2,  AIMAG( ckm ) * s2,  KIND = wp ) -                &
     804                        CMPLX( REAL( cjm ) * s72, AIMAG( cjm ) * s72, KIND = wp )
     805                   array(k2) = ck + CMPLX( -AIMAG( cj ), REAL( cj ), KIND = wp )
     806                   array(k3) = ck + CMPLX( AIMAG( cj ), - REAL( cj ), KIND = wp )
     807                   kk = k4 + kspan
     808                   IF ( kk >= nn )  EXIT
     809                END DO
     810                kk = kk - nn
     811                IF ( kk > kspan )  EXIT
     812             END DO
     813
     814          CASE default
     815             IF ( k /= jf )  THEN
     816                jf = k
     817                s1 = pi2 / k
     818                c1 = COS( s1 )
     819                s1 = SIN( s1 )
     820                cosine(jf) = 1.0_wp
     821                sine(jf) = 0.0_wp
     822                j = 1
     823                DO
     824                   cosine(j) = cosine(k) * c1 + sine(k) * s1
     825                   sine(j) = cosine(k) * s1 - sine(k) * c1
     826                   k = k - 1
     827                   cosine(k) = cosine(j)
     828                   sine(k) = - sine(j)
     829                   j = j + 1
     830                   IF ( j >= k )  EXIT
     831                END DO
     832             END IF
     833             DO
     834                DO
     835                   k1 = kk
     836                   k2 = kk + ispan
     837                   cc = array(kk)
     838                   ck = cc
     839                   j = 1
     840                   k1 = k1 + kspan
     841                   DO
     842                      k2 = k2 - kspan
     843                      j = j + 1
     844                      ctmp(j) = array(k1) + array(k2)
     845                      ck = ck + ctmp(j)
     846                      j = j + 1
     847                      ctmp(j) = array(k1) - array(k2)
     848                      k1 = k1 + kspan
     849                      IF ( k1 >= k2 )  EXIT
     850                   END DO
     851                   array(kk) = ck
     852                   k1 = kk
     853                   k2 = kk + ispan
     854                   j = 1
     855                   DO
     856                      k1 = k1 + kspan
     857                      k2 = k2 - kspan
     858                      jj = j
     859                      ck = cc
     860                      cj = ( 0.0_wp, 0.0_wp )
     861                      k = 1
     862                      DO
     863                         k = k + 1
     864                         ck = ck + CMPLX( REAL( ctmp(k) ) * cosine(jj), AIMAG( ctmp(k) ) *         &
     865                                          cosine(jj), KIND = wp )
     866                         k = k + 1
     867                         cj = cj + CMPLX( REAL( ctmp(k) ) * sine(jj), AIMAG( ctmp(k) ) * sine(jj), &
     868                                          KIND = wp )
     869                         jj = jj + j
     870                         IF ( jj > jf )  jj = jj - jf
     871                         IF ( k >= jf )  EXIT
     872                      END DO
     873                      k = jf - j
     874                      array(k1) = ck + CMPLX( - AIMAG( cj ), REAL( cj ), KIND = wp )
     875                      array(k2) = ck + CMPLX( AIMAG( cj ), -REAL( cj ), KIND = wp )
     876                      j = j + 1
     877                      IF ( j >= k )  EXIT
     878                   END DO
     879                   kk = kk + ispan
     880                   IF ( kk > nn )  EXIT
     881                END DO
     882                kk = kk - nn
     883                IF ( kk > kspan )  EXIT
     884             END DO
     885
     886          END SELECT
     887!
     888!--       Multiply by rotation factor (except for factors of 2 and 4)
     889          IF ( ii == nfactor )  RETURN
     890          kk = jc + 1
     891          DO
     892             c2 = 1.0_wp - cd
     893             s1 = sd
     894             DO
     895                c1 = c2
     896                s2 = s1
     897                kk = kk + kspan
     898                DO
     899                   DO
     900                      array(kk) = CMPLX( c2, s2, KIND = wp ) * array(kk)
     901                      kk = kk + ispan
     902                      IF ( kk > nt )  EXIT
     903                   END DO
     904                   ak = s1 * s2
     905                   s2 = s1 * c2 + c1 * s2
     906                   c2 = c1 * c2 - ak
     907                   kk = kk - nt + kspan
     908                   IF ( kk > ispan )  EXIT
     909                END DO
     910                c2 = c1 - ( cd * c1 + sd * s1 )
     911                s1 = s1 + sd * c1 - cd * s1
     912                c1 = 2.0_wp - ( c2 * c2 + s1 * s1 )
     913                s1 = s1 * c1
     914                c2 = c2 * c1
     915                kk = kk - ispan + jc
     916                IF ( kk > kspan )  EXIT
     917             END DO
     918             kk = kk - kspan + jc + 1
     919             IF ( kk > jc + jc )  EXIT
     920          END DO
     921
     922       END SELECT
     923    END DO
     924 END SUBROUTINE transform
     925
     926
     927!--------------------------------------------------------------------------------------------------!
    948928! Description:
    949929! ------------
    950930!> @todo Missing subroutine description.
    951 !------------------------------------------------------------------------------!
    952     SUBROUTINE permute(array, ntotal, npass, nspan, &
    953          factor, nfactor, nsquare, maxfactor, &
    954          ctmp, perm)
    955 !
    956 !--   Formal parameters
    957       COMPLEX(wp),  DIMENSION(*), INTENT(IN OUT):: array
    958       COMPLEX(wp),  DIMENSION(*), INTENT(OUT)   :: ctmp
    959       INTEGER(iwp),               INTENT(IN)    :: ntotal, npass, nspan
    960       INTEGER(iwp), DIMENSION(*), INTENT(IN OUT):: factor
    961       INTEGER(iwp),               INTENT(IN)    :: nfactor, nsquare
    962       INTEGER(iwp),               INTENT(IN)    :: maxfactor
    963       INTEGER(iwp), DIMENSION(*), INTENT(OUT)   :: perm
    964 !
    965 !--   Local scalars
    966       COMPLEX(wp) :: ck
    967       INTEGER(iwp):: ii, ispan
    968       INTEGER(iwp):: j, jc, jj
    969       INTEGER(iwp):: k, kk, kspan, kt, k1, k2, k3
    970       INTEGER(iwp):: nn, nt
    971 !
    972 !--   Permute the results to normal order---done in two stages
    973 !--   Permutation for square factors of n
    974 
    975       nt = ntotal
    976       nn = nt - 1
    977       kt = nsquare
    978       kspan = nspan
    979       jc = nspan / npass
    980 
    981       perm (1) = nspan
    982       IF (kt > 0) THEN
    983          k = kt + kt + 1
    984          IF (nfactor < k) k = k - 1
    985          j = 1
    986          perm (k + 1) = jc
    987          DO
    988             perm (j + 1) = perm (j) / factor(j)
    989             perm (k) = perm (k + 1) * factor(j)
    990             j = j + 1
    991             k = k - 1
    992             IF (j >= k) EXIT
    993          END DO
    994          k3 = perm (k + 1)
    995          kspan = perm (2)
    996          kk = jc + 1
    997          k2 = kspan + 1
    998          j = 1
    999 
    1000          IF (npass /= ntotal) THEN
    1001             permute_multi: DO
    1002                DO
    1003                   DO
    1004                      k = kk + jc
    1005                      DO
    1006 !
    1007 !--                     Swap array(kk) <> array(k2)
    1008                         ck = array(kk)
    1009                         array(kk) = array(k2)
    1010                         array(k2) = ck
    1011                         kk = kk + 1
    1012                         k2 = k2 + 1
    1013                         IF (kk >= k) EXIT
    1014                      END DO
    1015                      kk = kk + nspan - jc
    1016                      k2 = k2 + nspan - jc
    1017                      IF (kk >= nt) EXIT
    1018                   END DO
    1019                   kk = kk - nt + jc
    1020                   k2 = k2 - nt + kspan
    1021                   IF (k2 >= nspan) EXIT
    1022                END DO
    1023                DO
    1024                   DO
    1025                      k2 = k2 - perm (j)
    1026                      j = j + 1
    1027                      k2 = perm (j + 1) + k2
    1028                      IF (k2 <= perm (j)) EXIT
    1029                   END DO
    1030                   j = 1
    1031                   DO
    1032                      IF (kk < k2) CYCLE permute_multi
    1033                      kk = kk + jc
    1034                      k2 = k2 + kspan
    1035                      IF (k2 >= nspan) EXIT
    1036                   END DO
    1037                   IF (kk >= nspan) EXIT
    1038                END DO
    1039                EXIT
    1040             END DO permute_multi
    1041          ELSE
    1042             permute_single: DO
    1043                DO
    1044 !
    1045 !--               Swap array(kk) <> array(k2)
    1046                   ck = array(kk)
    1047                   array(kk) = array(k2)
    1048                   array(k2) = ck
    1049                   kk = kk + 1
    1050                   k2 = k2 + kspan
    1051                   IF (k2 >= nspan) EXIT
    1052                END DO
    1053                DO
    1054                   DO
    1055                      k2 = k2 - perm (j)
    1056                      j = j + 1
    1057                      k2 = perm (j + 1) + k2
    1058                      IF (k2 <= perm (j)) EXIT
    1059                   END DO
    1060                   j = 1
    1061                   DO
    1062                      IF (kk < k2) CYCLE permute_single
    1063                      kk = kk + 1
    1064                      k2 = k2 + kspan
    1065                      IF (k2 >= nspan) EXIT
    1066                   END DO
    1067                   IF (kk >= nspan) EXIT
    1068                END DO
    1069                EXIT
    1070             END DO permute_single
    1071          END IF
    1072          jc = k3
    1073       END IF
    1074 
    1075       IF (ISHFT(kt, 1) + 1 >= nfactor) RETURN
    1076 
    1077       ispan = perm (kt + 1)
    1078 !
    1079 !--   Permutation for square-free factors of n
    1080       j = nfactor - kt
    1081       factor(j + 1) = 1
    1082       DO
    1083          factor(j) = factor(j) * factor(j+1)
    1084          j = j - 1
    1085          IF (j == kt) EXIT
    1086       END DO
    1087       kt = kt + 1
    1088       nn = factor(kt) - 1
    1089       j = 0
    1090       jj = 0
    1091       DO
    1092          k = kt + 1
    1093          k2 = factor(kt)
    1094          kk = factor(k)
    1095          j = j + 1
    1096          IF (j > nn) EXIT !-- exit infinite loop
    1097          jj = jj + kk
    1098          DO WHILE (jj >= k2)
    1099             jj = jj - k2
    1100             k2 = kk
    1101             k = k + 1
    1102             kk = factor(k)
    1103             jj = jj + kk
    1104          END DO
    1105          perm (j) = jj
    1106       END DO
    1107 !
    1108 !--   Determine the permutation cycles of length greater than 1
    1109       j = 0
    1110       DO
    1111          DO
    1112             j = j + 1
    1113             kk = perm(j)
    1114             IF (kk >= 0) EXIT
    1115          END DO
    1116          IF (kk /= j) THEN
    1117             DO
    1118                k = kk
    1119                kk = perm (k)
    1120                perm (k) = -kk
    1121                IF (kk == j) EXIT
    1122             END DO
    1123             k3 = kk
    1124          ELSE
    1125             perm (j) = -j
    1126             IF (j == nn) EXIT !-- exit infinite loop
    1127          END IF
    1128       END DO
    1129 !
    1130 !--   Reorder a and b, following the permutation cycles
    1131       DO
    1132          j = k3 + 1
    1133          nt = nt - ispan
    1134          ii = nt - 1 + 1
    1135          IF (nt < 0) EXIT !-- exit infinite loop
    1136          DO
    1137             DO
    1138                j = j-1
    1139                IF (perm(j) >= 0) EXIT
    1140             END DO
    1141             jj = jc
    1142             DO
    1143                kspan = jj
    1144                IF (jj > maxfactor) kspan = maxfactor
    1145                jj = jj - kspan
    1146                k = perm(j)
    1147                kk = jc * k + ii + jj
    1148                k1 = kk + kspan
    1149                k2 = 0
    1150                DO
    1151                   k2 = k2 + 1
    1152                   ctmp(k2) = array(k1)
    1153                   k1 = k1 - 1
    1154                   IF (k1 == kk) EXIT
    1155                END DO
    1156                DO
    1157                   k1 = kk + kspan
    1158                   k2 = k1 - jc * (k + perm(k))
    1159                   k = -perm(k)
    1160                   DO
    1161                      array(k1) = array(k2)
    1162                      k1 = k1 - 1
    1163                      k2 = k2 - 1
    1164                      IF (k1 == kk) EXIT
    1165                   END DO
    1166                   kk = k2
    1167                   IF (k == j) EXIT
    1168                END DO
    1169                k1 = kk + kspan
    1170                k2 = 0
    1171                DO
    1172                   k2 = k2 + 1
    1173                   array(k1) = ctmp(k2)
    1174                   k1 = k1 - 1
    1175                   IF (k1 == kk) EXIT
    1176                END DO
    1177                IF (jj == 0) EXIT
    1178             END DO
    1179             IF (j == 1) EXIT
    1180          END DO
    1181       END DO
    1182 
    1183     END SUBROUTINE permute
     931!--------------------------------------------------------------------------------------------------!
     932 SUBROUTINE permute( array, ntotal, npass, nspan, factor, nfactor, nsquare, maxfactor, ctmp, perm )
     933!
     934!-- Formal parameters
     935    COMPLEX(wp),  DIMENSION(*), INTENT(IN OUT) ::  array                 !<
     936    COMPLEX(wp),  DIMENSION(*), INTENT(OUT)    ::  ctmp                  !<
     937    INTEGER(iwp),               INTENT(IN)     ::  ntotal, npass, nspan  !<
     938    INTEGER(iwp),               INTENT(IN)     ::  nfactor, nsquare      !<
     939    INTEGER(iwp),               INTENT(IN)     ::  maxfactor             !<
     940    INTEGER(iwp), DIMENSION(*), INTENT(IN OUT) ::  factor                !<
     941    INTEGER(iwp), DIMENSION(*), INTENT(OUT)    ::  perm                  !<
     942!
     943!-- Local scalars
     944    COMPLEX(wp)  ::  ck                            !<
     945    INTEGER(iwp) ::  ii, ispan                     !<
     946    INTEGER(iwp) ::  j, jc, jj                     !<
     947    INTEGER(iwp) ::  k, kk, kspan, kt, k1, k2, k3  !<
     948    INTEGER(iwp) ::  nn, nt                        !<
     949!
     950!-- Permute the results to normal order---done in two stages
     951!-- Permutation for square factors of n
     952
     953    nt = ntotal
     954    nn = nt - 1
     955    kt = nsquare
     956    kspan = nspan
     957    jc = nspan / npass
     958
     959    perm (1) = nspan
     960    IF ( kt > 0 )  THEN
     961       k = kt + kt + 1
     962       IF ( nfactor < k )  k = k - 1
     963       j = 1
     964       perm(k + 1) = jc
     965       DO
     966          perm(j + 1) = perm(j) / factor(j)
     967          perm(k) = perm(k + 1) * factor(j)
     968          j = j + 1
     969          k = k - 1
     970          IF ( j >= k )  EXIT
     971       END DO
     972       k3 = perm(k + 1)
     973       kspan = perm(2)
     974       kk = jc + 1
     975       k2 = kspan + 1
     976       j = 1
     977
     978       IF ( npass /= ntotal )  THEN
     979          permute_multi: DO
     980             DO
     981                DO
     982                   k = kk + jc
     983                   DO
     984!
     985!--                   Swap array(kk) <> array(k2)
     986                      ck = array(kk)
     987                      array(kk) = array(k2)
     988                      array(k2) = ck
     989                      kk = kk + 1
     990                      k2 = k2 + 1
     991                      IF ( kk >= k )  EXIT
     992                   END DO
     993                   kk = kk + nspan - jc
     994                   k2 = k2 + nspan - jc
     995                   IF ( kk >= nt )  EXIT
     996                END DO
     997                kk = kk - nt + jc
     998                k2 = k2 - nt + kspan
     999                IF ( k2 >= nspan )  EXIT
     1000             END DO
     1001             DO
     1002                DO
     1003                   k2 = k2 - perm(j)
     1004                   j = j + 1
     1005                   k2 = perm(j + 1) + k2
     1006                   IF ( k2 <= perm(j) )  EXIT
     1007                END DO
     1008                j = 1
     1009                DO
     1010                   IF ( kk < k2 )  CYCLE permute_multi
     1011                   kk = kk + jc
     1012                   k2 = k2 + kspan
     1013                   IF ( k2 >= nspan )  EXIT
     1014                END DO
     1015                IF ( kk >= nspan )  EXIT
     1016             END DO
     1017             EXIT
     1018          END DO permute_multi
     1019       ELSE
     1020          permute_single: DO
     1021             DO
     1022!
     1023!--             Swap array(kk) <> array(k2)
     1024                ck = array(kk)
     1025                array(kk) = array(k2)
     1026                array(k2) = ck
     1027                kk = kk + 1
     1028                k2 = k2 + kspan
     1029                IF ( k2 >= nspan )  EXIT
     1030             END DO
     1031             DO
     1032                DO
     1033                   k2 = k2 - perm(j)
     1034                   j = j + 1
     1035                   k2 = perm(j + 1) + k2
     1036                   IF ( k2 <= perm(j) )  EXIT
     1037                END DO
     1038                j = 1
     1039                DO
     1040                   IF ( kk < k2 )  CYCLE permute_single
     1041                   kk = kk + 1
     1042                   k2 = k2 + kspan
     1043                   IF ( k2 >= nspan )  EXIT
     1044                END DO
     1045                IF ( kk >= nspan )  EXIT
     1046             END DO
     1047             EXIT
     1048          END DO permute_single
     1049       END IF
     1050       jc = k3
     1051    END IF
     1052
     1053    IF ( ISHFT( kt, 1 ) + 1 >= nfactor )  RETURN
     1054
     1055    ispan = perm(kt + 1)
     1056!
     1057!-- Permutation for square-free factors of n
     1058    j = nfactor - kt
     1059    factor( j + 1 ) = 1
     1060    DO
     1061       factor(j) = factor(j) * factor(j+1)
     1062       j = j - 1
     1063       IF ( j == kt )  EXIT
     1064    END DO
     1065    kt = kt + 1
     1066    nn = factor( kt ) - 1
     1067    j = 0
     1068    jj = 0
     1069    DO
     1070       k = kt + 1
     1071       k2 = factor(kt)
     1072       kk = factor(k)
     1073       j = j + 1
     1074       IF ( j > nn )  EXIT !-- Exit infinite loop
     1075       jj = jj + kk
     1076       DO WHILE ( jj >= k2 )
     1077          jj = jj - k2
     1078          k2 = kk
     1079          k = k + 1
     1080          kk = factor(k)
     1081          jj = jj + kk
     1082       END DO
     1083       perm(j) = jj
     1084    END DO
     1085!
     1086!-- Determine the permutation cycles of length greater than 1
     1087    j = 0
     1088    DO
     1089       DO
     1090          j = j + 1
     1091          kk = perm(j)
     1092          IF ( kk >= 0 )  EXIT
     1093       END DO
     1094       IF ( kk /= j )  THEN
     1095          DO
     1096             k = kk
     1097             kk = perm(k)
     1098             perm(k) = - kk
     1099             IF ( kk == j )  EXIT
     1100          END DO
     1101          k3 = kk
     1102       ELSE
     1103          perm(j) = - j
     1104          IF ( j == nn )  EXIT !-- Exit infinite loop
     1105       END IF
     1106    END DO
     1107!
     1108!-- Reorder a and b, following the permutation cycles
     1109    DO
     1110       j = k3 + 1
     1111       nt = nt - ispan
     1112       ii = nt - 1 + 1
     1113       IF ( nt < 0 )  EXIT !-- Exit infinite loop
     1114       DO
     1115          DO
     1116             j = j - 1
     1117             IF ( perm(j) >= 0 )  EXIT
     1118          END DO
     1119          jj = jc
     1120          DO
     1121             kspan = jj
     1122             IF ( jj > maxfactor )  kspan = maxfactor
     1123             jj = jj - kspan
     1124             k = perm(j)
     1125             kk = jc * k + ii + jj
     1126             k1 = kk + kspan
     1127             k2 = 0
     1128             DO
     1129                k2 = k2 + 1
     1130                ctmp(k2) = array(k1)
     1131                k1 = k1 - 1
     1132                IF ( k1 == kk )  EXIT
     1133             END DO
     1134             DO
     1135                k1 = kk + kspan
     1136                k2 = k1 - jc * ( k + perm(k) )
     1137                k = - perm(k)
     1138                DO
     1139                   array(k1) = array(k2)
     1140                   k1 = k1 - 1
     1141                   k2 = k2 - 1
     1142                   IF ( k1 == kk )  EXIT
     1143                END DO
     1144                kk = k2
     1145                IF ( k == j )  EXIT
     1146             END DO
     1147             k1 = kk + kspan
     1148             k2 = 0
     1149             DO
     1150                k2 = k2 + 1
     1151                array(k1) = ctmp(k2)
     1152                k1 = k1 - 1
     1153                IF ( k1 == kk )  EXIT
     1154             END DO
     1155             IF ( jj == 0 )  EXIT
     1156          END DO
     1157          IF ( j == 1 )  EXIT
     1158       END DO
     1159    END DO
     1160
     1161 END SUBROUTINE permute
    11841162
    11851163 END SUBROUTINE fftradix
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |