|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 !Copyright by Stefan Goedecker, Cornell, Ithaca, USA, March 25, 1994 00003 !modified by Stefan Goedecker, Stuttgart, Germany, October 15, 1995 00004 !Commercial use is prohibited without the explicit permission of the author. 00005 !-----------------------------------------------------------------------------! 00006 00007 ! ***************************************************************************** 00008 SUBROUTINE fftstp ( mm, nfft, m, nn, n, zin, zout, trig, now, after, before, isign ) 00009 00010 USE fft_kinds, ONLY : dp 00011 INTEGER, INTENT ( IN ) :: mm, nfft, m, nn, n, now, after, before, isign 00012 REAL ( dp ), DIMENSION ( 2, 1024 ), INTENT ( IN ) :: trig 00013 REAL ( dp ), DIMENSION ( 2, mm, m ), INTENT ( IN ) :: zin 00014 REAL ( dp ), DIMENSION ( 2, nn, n ), INTENT ( INOUT ) :: zout 00015 00016 INTEGER :: atn, atb, ia, ib, nin1, nin2, nin3, nin4, nin5, nin6, nin7, nin8 00017 INTEGER :: nout1, nout2, nout3, nout4, nout5, nout6, nout7, nout8, 00018 j, ias, itt, itrig 00019 REAL ( dp ) :: s, s1, s2, s3, s4, s5, s6, s7, s8, 00020 r, r1, r2, r3, r4, r5, r6, r7, r8, cr2, cr3, cr4, cr5, 00021 ci2, ci3, ci4, ci5, ur1, ur2, ur3, ui1, ui2, ui3, 00022 vr1, vr2, vr3, vi1, vi2, vi3, cm, cp, dm, dbl, 00023 am, ap, bm, bp, bbs, s25, s34, r34, r25, sin2, sin4 00024 REAL ( dp ), PARAMETER :: rt2i = 0.7071067811865475_dp ! sqrt(0.5) 00025 REAL ( dp ), PARAMETER :: bb = 0.8660254037844387_dp ! sqrt(3)/2 00026 REAL ( dp ), PARAMETER :: cos2 = 0.3090169943749474_dp ! cos(2*pi/5) 00027 REAL ( dp ), PARAMETER :: cos4 = - 0.8090169943749474_dp ! cos(4*pi/5) 00028 REAL ( dp ), PARAMETER :: sin2p = 0.9510565162951536_dp ! sin(2*pi/5) 00029 REAL ( dp ), PARAMETER :: sin4p = 0.5877852522924731_dp ! sin(4*pi/5) 00030 00031 !-----------------------------------------------------------------------------! 00032 00033 atn = after * now 00034 atb = after * before 00035 00036 IF ( now == 4 ) THEN 00037 IF ( isign == 1 ) THEN 00038 ia = 1 00039 nin1 = ia - after 00040 nout1 = ia - atn 00041 DO ib = 1, before 00042 nin1 = nin1 + after 00043 nin2 = nin1 + atb 00044 nin3 = nin2 + atb 00045 nin4 = nin3 + atb 00046 nout1 = nout1 + atn 00047 nout2 = nout1 + after 00048 nout3 = nout2 + after 00049 nout4 = nout3 + after 00050 DO j = 1, nfft 00051 r1 = zin ( 1, j, nin1 ) 00052 s1 = zin ( 2, j, nin1 ) 00053 r2 = zin ( 1, j, nin2 ) 00054 s2 = zin ( 2, j, nin2 ) 00055 r3 = zin ( 1, j, nin3 ) 00056 s3 = zin ( 2, j, nin3 ) 00057 r4 = zin ( 1, j, nin4 ) 00058 s4 = zin ( 2, j, nin4 ) 00059 r = r1 + r3 00060 s = r2 + r4 00061 zout ( 1, j, nout1 ) = r + s 00062 zout ( 1, j, nout3 ) = r - s 00063 r = r1 - r3 00064 s = s2 - s4 00065 zout ( 1, j, nout2 ) = r - s 00066 zout ( 1, j, nout4 ) = r + s 00067 r = s1 + s3 00068 s = s2 + s4 00069 zout ( 2, j, nout1 ) = r + s 00070 zout ( 2, j, nout3 ) = r - s 00071 r = s1 - s3 00072 s = r2 - r4 00073 zout ( 2, j, nout2 ) = r + s 00074 zout ( 2, j, nout4 ) = r - s 00075 END DO 00076 END DO 00077 DO ia = 2, after 00078 ias = ia - 1 00079 IF ( 2*ias == after ) THEN 00080 nin1 = ia - after 00081 nout1 = ia - atn 00082 DO ib = 1, before 00083 nin1 = nin1 + after 00084 nin2 = nin1 + atb 00085 nin3 = nin2 + atb 00086 nin4 = nin3 + atb 00087 nout1 = nout1 + atn 00088 nout2 = nout1 + after 00089 nout3 = nout2 + after 00090 nout4 = nout3 + after 00091 DO j = 1, nfft 00092 r1 = zin ( 1, j, nin1 ) 00093 s1 = zin ( 2, j, nin1 ) 00094 r = zin ( 1, j, nin2 ) 00095 s = zin ( 2, j, nin2 ) 00096 r2 = ( r - s ) * rt2i 00097 s2 = ( r + s ) * rt2i 00098 r3 = - zin ( 2, j, nin3 ) 00099 s3 = zin ( 1, j, nin3 ) 00100 r = zin ( 1, j, nin4 ) 00101 s = zin ( 2, j, nin4 ) 00102 r4 = - ( r + s ) * rt2i 00103 s4 = ( r - s ) * rt2i 00104 r = r1 + r3 00105 s = r2 + r4 00106 zout ( 1, j, nout1 ) = r + s 00107 zout ( 1, j, nout3 ) = r - s 00108 r = r1 - r3 00109 s = s2 - s4 00110 zout ( 1, j, nout2 ) = r - s 00111 zout ( 1, j, nout4 ) = r + s 00112 r = s1 + s3 00113 s = s2 + s4 00114 zout ( 2, j, nout1 ) = r + s 00115 zout ( 2, j, nout3 ) = r - s 00116 r = s1 - s3 00117 s = r2 - r4 00118 zout ( 2, j, nout2 ) = r + s 00119 zout ( 2, j, nout4 ) = r - s 00120 END DO 00121 END DO 00122 ELSE 00123 itt = ias * before 00124 itrig = itt + 1 00125 cr2 = trig ( 1, itrig ) 00126 ci2 = trig ( 2, itrig ) 00127 itrig = itrig + itt 00128 cr3 = trig ( 1, itrig ) 00129 ci3 = trig ( 2, itrig ) 00130 itrig = itrig + itt 00131 cr4 = trig ( 1, itrig ) 00132 ci4 = trig ( 2, itrig ) 00133 nin1 = ia - after 00134 nout1 = ia - atn 00135 DO ib = 1, before 00136 nin1 = nin1 + after 00137 nin2 = nin1 + atb 00138 nin3 = nin2 + atb 00139 nin4 = nin3 + atb 00140 nout1 = nout1 + atn 00141 nout2 = nout1 + after 00142 nout3 = nout2 + after 00143 nout4 = nout3 + after 00144 DO j = 1, nfft 00145 r1 = zin ( 1, j, nin1 ) 00146 s1 = zin ( 2, j, nin1 ) 00147 r = zin ( 1, j, nin2 ) 00148 s = zin ( 2, j, nin2 ) 00149 r2 = r * cr2 - s * ci2 00150 s2 = r * ci2 + s * cr2 00151 r = zin ( 1, j, nin3 ) 00152 s = zin ( 2, j, nin3 ) 00153 r3 = r * cr3 - s * ci3 00154 s3 = r * ci3 + s * cr3 00155 r = zin ( 1, j, nin4 ) 00156 s = zin ( 2, j, nin4 ) 00157 r4 = r * cr4 - s * ci4 00158 s4 = r * ci4 + s * cr4 00159 r = r1 + r3 00160 s = r2 + r4 00161 zout ( 1, j, nout1 ) = r + s 00162 zout ( 1, j, nout3 ) = r - s 00163 r = r1 - r3 00164 s = s2 - s4 00165 zout ( 1, j, nout2 ) = r - s 00166 zout ( 1, j, nout4 ) = r + s 00167 r = s1 + s3 00168 s = s2 + s4 00169 zout ( 2, j, nout1 ) = r + s 00170 zout ( 2, j, nout3 ) = r - s 00171 r = s1 - s3 00172 s = r2 - r4 00173 zout ( 2, j, nout2 ) = r + s 00174 zout ( 2, j, nout4 ) = r - s 00175 END DO 00176 END DO 00177 END IF 00178 END DO 00179 ELSE 00180 ia = 1 00181 nin1 = ia - after 00182 nout1 = ia - atn 00183 DO ib = 1, before 00184 nin1 = nin1 + after 00185 nin2 = nin1 + atb 00186 nin3 = nin2 + atb 00187 nin4 = nin3 + atb 00188 nout1 = nout1 + atn 00189 nout2 = nout1 + after 00190 nout3 = nout2 + after 00191 nout4 = nout3 + after 00192 DO j = 1, nfft 00193 r1 = zin ( 1, j, nin1 ) 00194 s1 = zin ( 2, j, nin1 ) 00195 r2 = zin ( 1, j, nin2 ) 00196 s2 = zin ( 2, j, nin2 ) 00197 r3 = zin ( 1, j, nin3 ) 00198 s3 = zin ( 2, j, nin3 ) 00199 r4 = zin ( 1, j, nin4 ) 00200 s4 = zin ( 2, j, nin4 ) 00201 r = r1 + r3 00202 s = r2 + r4 00203 zout ( 1, j, nout1 ) = r + s 00204 zout ( 1, j, nout3 ) = r - s 00205 r = r1 - r3 00206 s = s2 - s4 00207 zout ( 1, j, nout2 ) = r + s 00208 zout ( 1, j, nout4 ) = r - s 00209 r = s1 + s3 00210 s = s2 + s4 00211 zout ( 2, j, nout1 ) = r + s 00212 zout ( 2, j, nout3 ) = r - s 00213 r = s1 - s3 00214 s = r2 - r4 00215 zout ( 2, j, nout2 ) = r - s 00216 zout ( 2, j, nout4 ) = r + s 00217 END DO 00218 END DO 00219 DO ia = 2, after 00220 ias = ia - 1 00221 IF ( 2 * ias == after ) THEN 00222 nin1 = ia - after 00223 nout1 = ia - atn 00224 DO ib = 1, before 00225 nin1 = nin1 + after 00226 nin2 = nin1 + atb 00227 nin3 = nin2 + atb 00228 nin4 = nin3 + atb 00229 nout1 = nout1 + atn 00230 nout2 = nout1 + after 00231 nout3 = nout2 + after 00232 nout4 = nout3 + after 00233 DO j = 1, nfft 00234 r1 = zin ( 1, j, nin1 ) 00235 s1 = zin ( 2, j, nin1 ) 00236 r = zin ( 1, j, nin2 ) 00237 s = zin ( 2, j, nin2 ) 00238 r2 = ( r + s ) * rt2i 00239 s2 = ( s - r ) * rt2i 00240 r3 = zin ( 2, j, nin3 ) 00241 s3 = - zin ( 1, j, nin3 ) 00242 r = zin ( 1, j, nin4 ) 00243 s = zin ( 2, j, nin4 ) 00244 r4 = ( s - r ) * rt2i 00245 s4 = - ( r + s ) * rt2i 00246 r = r1 + r3 00247 s = r2 + r4 00248 zout ( 1, j, nout1 ) = r + s 00249 zout ( 1, j, nout3 ) = r - s 00250 r = r1 - r3 00251 s = s2 - s4 00252 zout ( 1, j, nout2 ) = r + s 00253 zout ( 1, j, nout4 ) = r - s 00254 r =s1 + s3 00255 s =s2 + s4 00256 zout ( 2, j, nout1 ) = r + s 00257 zout ( 2, j, nout3 ) = r - s 00258 r = s1 - s3 00259 s = r2 - r4 00260 zout ( 2, j, nout2 ) = r - s 00261 zout ( 2, j, nout4 ) = r + s 00262 END DO 00263 END DO 00264 ELSE 00265 itt = ias * before 00266 itrig = itt + 1 00267 cr2 = trig ( 1, itrig ) 00268 ci2 = trig ( 2, itrig ) 00269 itrig = itrig + itt 00270 cr3 = trig ( 1, itrig ) 00271 ci3 = trig ( 2, itrig ) 00272 itrig = itrig + itt 00273 cr4 = trig ( 1, itrig ) 00274 ci4 = trig ( 2, itrig ) 00275 nin1 = ia - after 00276 nout1 = ia - atn 00277 DO ib = 1, before 00278 nin1 = nin1 + after 00279 nin2 = nin1 + atb 00280 nin3 = nin2 + atb 00281 nin4 = nin3 + atb 00282 nout1 = nout1 + atn 00283 nout2 = nout1 + after 00284 nout3 = nout2 + after 00285 nout4 = nout3 + after 00286 DO j = 1, nfft 00287 r1 = zin ( 1, j, nin1 ) 00288 s1 = zin ( 2, j, nin1 ) 00289 r = zin ( 1, j, nin2 ) 00290 s = zin ( 2, j, nin2 ) 00291 r2 = r * cr2 - s * ci2 00292 s2 = r * ci2 + s * cr2 00293 r = zin ( 1, j, nin3 ) 00294 s = zin ( 2, j, nin3 ) 00295 r3 = r * cr3 - s * ci3 00296 s3 = r * ci3 + s * cr3 00297 r = zin ( 1, j, nin4 ) 00298 s = zin ( 2, j, nin4 ) 00299 r4 = r * cr4 - s * ci4 00300 s4 = r * ci4 + s * cr4 00301 r = r1 + r3 00302 s = r2 + r4 00303 zout ( 1, j, nout1 ) = r + s 00304 zout ( 1, j, nout3 ) = r - s 00305 r = r1 - r3 00306 s = s2 - s4 00307 zout ( 1, j, nout2 ) = r + s 00308 zout ( 1, j, nout4 ) = r - s 00309 r = s1 + s3 00310 s = s2 + s4 00311 zout ( 2, j, nout1 ) = r + s 00312 zout ( 2, j, nout3 ) = r - s 00313 r = s1 - s3 00314 s = r2 - r4 00315 zout ( 2, j, nout2 ) = r - s 00316 zout ( 2, j, nout4 ) = r + s 00317 END DO 00318 END DO 00319 END IF 00320 END DO 00321 END IF 00322 ELSE IF ( now == 8 ) THEN 00323 IF ( isign == -1 ) THEN 00324 ia = 1 00325 nin1 = ia - after 00326 nout1 = ia - atn 00327 DO ib = 1, before 00328 nin1 = nin1 + after 00329 nin2 = nin1 + atb 00330 nin3 = nin2 + atb 00331 nin4 = nin3 + atb 00332 nin5 = nin4 + atb 00333 nin6 = nin5 + atb 00334 nin7 = nin6 + atb 00335 nin8 = nin7 + atb 00336 nout1 = nout1 + atn 00337 nout2 = nout1 + after 00338 nout3 = nout2 + after 00339 nout4 = nout3 + after 00340 nout5 = nout4 + after 00341 nout6 = nout5 + after 00342 nout7 = nout6 + after 00343 nout8 = nout7 + after 00344 DO j = 1, nfft 00345 r1 = zin ( 1, j, nin1 ) 00346 s1 = zin ( 2, j, nin1 ) 00347 r2 = zin ( 1, j, nin2 ) 00348 s2 = zin ( 2, j, nin2 ) 00349 r3 = zin ( 1, j, nin3 ) 00350 s3 = zin ( 2, j, nin3 ) 00351 r4 = zin ( 1, j, nin4 ) 00352 s4 = zin ( 2, j, nin4 ) 00353 r5 = zin ( 1, j, nin5 ) 00354 s5 = zin ( 2, j, nin5 ) 00355 r6 = zin ( 1, j, nin6 ) 00356 s6 = zin ( 2, j, nin6 ) 00357 r7 = zin ( 1, j, nin7 ) 00358 s7 = zin ( 2, j, nin7 ) 00359 r8 = zin ( 1, j, nin8 ) 00360 s8 = zin ( 2, j, nin8 ) 00361 r = r1 + r5 00362 s = r3 + r7 00363 ap = r + s 00364 am = r - s 00365 r = r2 + r6 00366 s = r4 + r8 00367 bp = r + s 00368 bm = r - s 00369 r = s1 + s5 00370 s = s3 + s7 00371 cp = r + s 00372 cm = r - s 00373 r = s2 + s6 00374 s = s4 + s8 00375 dbl = r + s 00376 dm = r - s 00377 zout ( 1, j, nout1 ) = ap + bp 00378 zout ( 2, j, nout1 ) = cp + dbl 00379 zout ( 1, j, nout5 ) = ap - bp 00380 zout ( 2, j, nout5 ) = cp - dbl 00381 zout ( 1, j, nout3 ) = am + dm 00382 zout ( 2, j, nout3 ) = cm - bm 00383 zout ( 1, j, nout7 ) = am - dm 00384 zout ( 2, j, nout7 ) = cm + bm 00385 r = r1 - r5 00386 s = s3 - s7 00387 ap = r + s 00388 am = r - s 00389 r = s1 - s5 00390 s = r3 - r7 00391 bp = r + s 00392 bm = r - s 00393 r = s4 - s8 00394 s = r2 - r6 00395 cp = r + s 00396 cm = r - s 00397 r = s2 - s6 00398 s = r4 - r8 00399 dbl = r + s 00400 dm = r - s 00401 r = ( cp + dm ) * rt2i 00402 s = (-cp + dm ) * rt2i 00403 cp = ( cm + dbl ) * rt2i 00404 dbl = ( cm - dbl ) * rt2i 00405 zout ( 1, j, nout2 ) = ap + r 00406 zout ( 2, j, nout2 ) = bm + s 00407 zout ( 1, j, nout6 ) = ap - r 00408 zout ( 2, j, nout6 ) = bm - s 00409 zout ( 1, j, nout4 ) = am + cp 00410 zout ( 2, j, nout4 ) = bp + dbl 00411 zout ( 1, j, nout8 ) = am - cp 00412 zout ( 2, j, nout8 ) = bp - dbl 00413 END DO 00414 END DO 00415 ELSE 00416 ia = 1 00417 nin1 = ia - after 00418 nout1 = ia - atn 00419 DO ib = 1, before 00420 nin1 = nin1 + after 00421 nin2 = nin1 + atb 00422 nin3 = nin2 + atb 00423 nin4 = nin3 + atb 00424 nin5 = nin4 + atb 00425 nin6 = nin5 + atb 00426 nin7 = nin6 + atb 00427 nin8 = nin7 + atb 00428 nout1 = nout1 + atn 00429 nout2 = nout1 + after 00430 nout3 = nout2 + after 00431 nout4 = nout3 + after 00432 nout5 = nout4 + after 00433 nout6 = nout5 + after 00434 nout7 = nout6 + after 00435 nout8 = nout7 + after 00436 DO j = 1, nfft 00437 r1 = zin ( 1, j, nin1 ) 00438 s1 = zin ( 2, j, nin1 ) 00439 r2 = zin ( 1, j, nin2 ) 00440 s2 = zin ( 2, j, nin2 ) 00441 r3 = zin ( 1, j, nin3 ) 00442 s3 = zin ( 2, j, nin3 ) 00443 r4 = zin ( 1, j, nin4 ) 00444 s4 = zin ( 2, j, nin4 ) 00445 r5 = zin ( 1, j, nin5 ) 00446 s5 = zin ( 2, j, nin5 ) 00447 r6 = zin ( 1, j, nin6 ) 00448 s6 = zin ( 2, j, nin6 ) 00449 r7 = zin ( 1, j, nin7 ) 00450 s7 = zin ( 2, j, nin7 ) 00451 r8 = zin ( 1, j, nin8 ) 00452 s8 = zin ( 2, j, nin8 ) 00453 r = r1 + r5 00454 s = r3 + r7 00455 ap = r + s 00456 am = r - s 00457 r = r2 + r6 00458 s = r4 + r8 00459 bp = r + s 00460 bm = r - s 00461 r = s1 + s5 00462 s = s3 + s7 00463 cp = r + s 00464 cm = r - s 00465 r = s2 + s6 00466 s = s4 + s8 00467 dbl = r + s 00468 dm = r - s 00469 zout ( 1, j, nout1 ) = ap + bp 00470 zout ( 2, j, nout1 ) = cp + dbl 00471 zout ( 1, j, nout5 ) = ap - bp 00472 zout ( 2, j, nout5 ) = cp - dbl 00473 zout ( 1, j, nout3 ) = am - dm 00474 zout ( 2, j, nout3 ) = cm + bm 00475 zout ( 1, j, nout7 ) = am + dm 00476 zout ( 2, j, nout7 ) = cm - bm 00477 r = r1 - r5 00478 s = -s3 + s7 00479 ap = r + s 00480 am = r - s 00481 r = s1 - s5 00482 s = r7 - r3 00483 bp = r + s 00484 bm = r - s 00485 r = -s4 + s8 00486 s = r2 - r6 00487 cp = r + s 00488 cm = r - s 00489 r = -s2 + s6 00490 s = r4 - r8 00491 dbl = r + s 00492 dm = r - s 00493 r = ( cp + dm ) * rt2i 00494 s = ( cp - dm ) * rt2i 00495 cp = ( cm + dbl ) * rt2i 00496 dbl = (-cm + dbl ) * rt2i 00497 zout ( 1, j, nout2 ) = ap + r 00498 zout ( 2, j, nout2 ) = bm + s 00499 zout ( 1, j, nout6 ) = ap - r 00500 zout ( 2, j, nout6 ) = bm - s 00501 zout ( 1, j, nout4 ) = am + cp 00502 zout ( 2, j, nout4 ) = bp + dbl 00503 zout ( 1, j, nout8 ) = am - cp 00504 zout ( 2, j, nout8 ) = bp - dbl 00505 END DO 00506 END DO 00507 END IF 00508 ELSE IF ( now == 3 ) THEN 00509 bbs = isign * bb 00510 ia = 1 00511 nin1 = ia - after 00512 nout1 = ia - atn 00513 DO ib = 1, before 00514 nin1 = nin1 + after 00515 nin2 = nin1 + atb 00516 nin3 = nin2 + atb 00517 nout1 = nout1 + atn 00518 nout2 = nout1 + after 00519 nout3 = nout2 + after 00520 DO j = 1, nfft 00521 r1 = zin ( 1, j, nin1 ) 00522 s1 = zin ( 2, j, nin1 ) 00523 r2 = zin ( 1, j, nin2 ) 00524 s2 = zin ( 2, j, nin2 ) 00525 r3 = zin ( 1, j, nin3 ) 00526 s3 = zin ( 2, j, nin3 ) 00527 r = r2 + r3 00528 s = s2 + s3 00529 zout ( 1, j, nout1 ) = r + r1 00530 zout ( 2, j, nout1 ) = s + s1 00531 r1 = r1 - 0.5_dp * r 00532 s1 = s1 - 0.5_dp * s 00533 r2 = bbs * ( r2 - r3 ) 00534 s2 = bbs * ( s2 - s3 ) 00535 zout ( 1, j, nout2 ) = r1 - s2 00536 zout ( 2, j, nout2 ) = s1 + r2 00537 zout ( 1, j, nout3 ) = r1 + s2 00538 zout ( 2, j, nout3 ) = s1 - r2 00539 END DO 00540 END DO 00541 DO ia = 2, after 00542 ias = ia - 1 00543 IF ( 4*ias == 3*after ) THEN 00544 IF ( isign == 1 ) THEN 00545 nin1 = ia - after 00546 nout1 = ia - atn 00547 DO ib = 1, before 00548 nin1 = nin1 + after 00549 nin2=nin1+atb 00550 nin3=nin2+atb 00551 nout1=nout1+atn 00552 nout2=nout1+after 00553 nout3=nout2+after 00554 DO j = 1, nfft 00555 r1 = zin ( 1, j, nin1 ) 00556 s1 = zin ( 2, j, nin1 ) 00557 r2 = -zin ( 2, j, nin2 ) 00558 s2 = zin ( 1, j, nin2 ) 00559 r3 = -zin ( 1, j, nin3 ) 00560 s3 = -zin ( 2, j, nin3 ) 00561 r = r2 + r3 00562 s = s2 + s3 00563 zout ( 1, j, nout1 ) = r + r1 00564 zout ( 2, j, nout1 ) = s + s1 00565 r1 = r1 - 0.5_dp * r 00566 s1 = s1 - 0.5_dp * s 00567 r2 = bbs*(r2-r3) 00568 s2 = bbs*(s2-s3) 00569 zout ( 1, j, nout2 ) = r1 - s2 00570 zout ( 2, j, nout2 ) = s1 + r2 00571 zout ( 1, j, nout3 ) = r1 + s2 00572 zout ( 2, j, nout3 ) = s1 - r2 00573 END DO 00574 END DO 00575 ELSE 00576 nin1 = ia - after 00577 nout1 = ia - atn 00578 DO ib = 1, before 00579 nin1 = nin1 + after 00580 nin2 = nin1 + atb 00581 nin3 = nin2 + atb 00582 nout1 = nout1 + atn 00583 nout2 = nout1 + after 00584 nout3 = nout2 + after 00585 DO j = 1, nfft 00586 r1 = zin ( 1, j, nin1 ) 00587 s1 = zin ( 2, j, nin1 ) 00588 r2 = zin ( 2, j, nin2 ) 00589 s2 = -zin ( 1, j, nin2 ) 00590 r3 = -zin ( 1, j, nin3 ) 00591 s3 = -zin ( 2, j, nin3 ) 00592 r = r2 + r3 00593 s = s2 + s3 00594 zout ( 1, j, nout1 ) = r + r1 00595 zout ( 2, j, nout1 ) = s + s1 00596 r1 = r1 - 0.5_dp * r 00597 s1 = s1 - 0.5_dp * s 00598 r2 = bbs * ( r2 - r3 ) 00599 s2 = bbs * ( s2 - s3 ) 00600 zout ( 1, j, nout2 ) = r1 - s2 00601 zout ( 2, j, nout2 ) = s1 + r2 00602 zout ( 1, j, nout3 ) = r1 + s2 00603 zout ( 2, j, nout3 ) = s1 - r2 00604 END DO 00605 END DO 00606 END IF 00607 ELSE IF ( 8 * ias == 3 * after ) THEN 00608 IF ( isign == 1 ) THEN 00609 nin1 = ia - after 00610 nout1 = ia - atn 00611 DO ib = 1, before 00612 nin1 = nin1 + after 00613 nin2 = nin1 + atb 00614 nin3 = nin2 + atb 00615 nout1 = nout1 + atn 00616 nout2 = nout1 + after 00617 nout3 = nout2 + after 00618 DO j = 1, nfft 00619 r1 = zin ( 1, j, nin1 ) 00620 s1 = zin ( 2, j, nin1 ) 00621 r = zin ( 1, j, nin2 ) 00622 s = zin ( 2, j, nin2 ) 00623 r2 = ( r - s ) * rt2i 00624 s2 = ( r + s ) * rt2i 00625 r3 = -zin ( 2, j, nin3 ) 00626 s3 = zin ( 1, j, nin3 ) 00627 r = r2 + r3 00628 s = s2 + s3 00629 zout ( 1, j, nout1 ) = r + r1 00630 zout ( 2, j, nout1 ) = s + s1 00631 r1 = r1 - 0.5_dp * r 00632 s1 = s1 - 0.5_dp * s 00633 r2 = bbs * ( r2 - r3 ) 00634 s2 = bbs * ( s2 - s3 ) 00635 zout ( 1, j, nout2 ) = r1 - s2 00636 zout ( 2, j, nout2 ) = s1 + r2 00637 zout ( 1, j, nout3 ) = r1 + s2 00638 zout ( 2, j, nout3 ) = s1 - r2 00639 END DO 00640 END DO 00641 ELSE 00642 nin1 = ia - after 00643 nout1 = ia - atn 00644 DO ib = 1, before 00645 nin1 = nin1 + after 00646 nin2 = nin1 + atb 00647 nin3 = nin2 + atb 00648 nout1 = nout1 + atn 00649 nout2 = nout1 + after 00650 nout3 = nout2 + after 00651 DO j = 1, nfft 00652 r1 = zin ( 1, j, nin1 ) 00653 s1 = zin ( 2, j, nin1 ) 00654 r = zin ( 1, j, nin2 ) 00655 s = zin ( 2, j, nin2 ) 00656 r2 = ( r + s ) * rt2i 00657 s2 = ( -r + s ) * rt2i 00658 r3 = zin ( 2, j, nin3 ) 00659 s3 = -zin ( 1, j, nin3 ) 00660 r = r2 + r3 00661 s = s2 + s3 00662 zout ( 1, j, nout1 ) = r + r1 00663 zout ( 2, j, nout1 ) = s + s1 00664 r1 = r1 - 0.5_dp * r 00665 s1 = s1 - 0.5_dp * s 00666 r2 = bbs * ( r2 - r3 ) 00667 s2 = bbs * ( s2 - s3 ) 00668 zout ( 1, j, nout2 ) = r1 - s2 00669 zout ( 2, j, nout2 ) = s1 + r2 00670 zout ( 1, j, nout3 ) = r1 + s2 00671 zout ( 2, j, nout3 ) = s1 - r2 00672 END DO 00673 END DO 00674 END IF 00675 ELSE 00676 itt = ias * before 00677 itrig = itt + 1 00678 cr2 = trig ( 1, itrig ) 00679 ci2 = trig ( 2, itrig ) 00680 itrig = itrig + itt 00681 cr3 = trig ( 1, itrig ) 00682 ci3 = trig ( 2, itrig ) 00683 nin1 = ia - after 00684 nout1 = ia - atn 00685 DO ib = 1, before 00686 nin1 = nin1 + after 00687 nin2 = nin1 + atb 00688 nin3 = nin2 + atb 00689 nout1 = nout1 + atn 00690 nout2 = nout1 + after 00691 nout3 = nout2 + after 00692 DO j = 1, nfft 00693 r1 = zin ( 1, j, nin1 ) 00694 s1 = zin ( 2, j, nin1 ) 00695 r = zin ( 1, j, nin2 ) 00696 s = zin ( 2, j, nin2 ) 00697 r2 = r * cr2 - s * ci2 00698 s2 = r * ci2 + s * cr2 00699 r = zin ( 1, j, nin3 ) 00700 s = zin ( 2, j, nin3 ) 00701 r3 = r * cr3 - s * ci3 00702 s3 = r * ci3 + s * cr3 00703 r = r2 + r3 00704 s = s2 + s3 00705 zout ( 1, j, nout1 ) = r + r1 00706 zout ( 2, j, nout1 ) = s + s1 00707 r1 = r1 - 0.5_dp * r 00708 s1 = s1 - 0.5_dp * s 00709 r2 = bbs * ( r2 - r3 ) 00710 s2 = bbs * ( s2 - s3 ) 00711 zout ( 1, j, nout2 ) = r1 - s2 00712 zout ( 2, j, nout2 ) = s1 + r2 00713 zout ( 1, j, nout3 ) = r1 + s2 00714 zout ( 2, j, nout3 ) = s1 - r2 00715 END DO 00716 END DO 00717 END IF 00718 END DO 00719 ELSE IF ( now == 5 ) THEN 00720 sin2 = isign * sin2p 00721 sin4 = isign * sin4p 00722 ia = 1 00723 nin1 = ia - after 00724 nout1 = ia - atn 00725 DO ib = 1, before 00726 nin1 = nin1 + after 00727 nin2 = nin1 + atb 00728 nin3 = nin2 + atb 00729 nin4 = nin3 + atb 00730 nin5 = nin4 + atb 00731 nout1 = nout1 + atn 00732 nout2 = nout1 + after 00733 nout3 = nout2 + after 00734 nout4 = nout3 + after 00735 nout5 = nout4 + after 00736 DO j = 1, nfft 00737 r1 = zin ( 1, j, nin1 ) 00738 s1 = zin ( 2, j, nin1 ) 00739 r2 = zin ( 1, j, nin2 ) 00740 s2 = zin ( 2, j, nin2 ) 00741 r3 = zin ( 1, j, nin3 ) 00742 s3 = zin ( 2, j, nin3 ) 00743 r4 = zin ( 1, j, nin4 ) 00744 s4 = zin ( 2, j, nin4 ) 00745 r5 = zin ( 1, j, nin5 ) 00746 s5 = zin ( 2, j, nin5 ) 00747 r25 = r2 + r5 00748 r34 = r3 + r4 00749 s25 = s2 - s5 00750 s34 = s3 - s4 00751 zout ( 1, j, nout1 ) = r1 + r25 + r34 00752 r = cos2 * r25 + cos4 * r34 + r1 00753 s = sin2 * s25 + sin4 * s34 00754 zout ( 1, j, nout2 ) = r - s 00755 zout ( 1, j, nout5 ) = r + s 00756 r = cos4 * r25 + cos2 * r34 + r1 00757 s = sin4 * s25 - sin2 * s34 00758 zout ( 1, j, nout3 ) = r - s 00759 zout ( 1, j, nout4 ) = r + s 00760 r25 = r2 - r5 00761 r34 = r3 - r4 00762 s25 = s2 + s5 00763 s34 = s3 + s4 00764 zout ( 2, j, nout1 ) = s1 + s25 + s34 00765 r = cos2 * s25 + cos4 * s34 + s1 00766 s = sin2 * r25 + sin4 * r34 00767 zout ( 2, j, nout2 ) = r + s 00768 zout ( 2, j, nout5 ) = r - s 00769 r = cos4 * s25 + cos2 * s34 + s1 00770 s = sin4 * r25 - sin2 * r34 00771 zout ( 2, j, nout3 ) = r + s 00772 zout ( 2, j, nout4 ) = r - s 00773 END DO 00774 END DO 00775 DO ia = 2, after 00776 ias = ia - 1 00777 IF ( 8 * ias == 5 * after ) THEN 00778 IF ( isign == 1 ) THEN 00779 nin1 = ia - after 00780 nout1 = ia - atn 00781 DO ib = 1, before 00782 nin1 = nin1 + after 00783 nin2 = nin1 + atb 00784 nin3 = nin2 + atb 00785 nin4 = nin3 + atb 00786 nin5 = nin4 + atb 00787 nout1 = nout1 + atn 00788 nout2 = nout1 + after 00789 nout3 = nout2 + after 00790 nout4 = nout3 + after 00791 nout5 = nout4 + after 00792 DO j = 1, nfft 00793 r1 = zin ( 1, j, nin1 ) 00794 s1 = zin ( 2, j, nin1 ) 00795 r = zin ( 1, j, nin2 ) 00796 s = zin ( 2, j, nin2 ) 00797 r2 = ( r - s ) * rt2i 00798 s2 = ( r + s ) * rt2i 00799 r3 = -zin ( 2, j, nin3 ) 00800 s3 = zin ( 1, j, nin3 ) 00801 r = zin ( 1, j, nin4 ) 00802 s = zin ( 2, j, nin4 ) 00803 r4 = -( r + s ) * rt2i 00804 s4 = ( r - s ) * rt2i 00805 r5 = -zin ( 1, j, nin5 ) 00806 s5 = -zin ( 2, j, nin5 ) 00807 r25 = r2 + r5 00808 r34 = r3 + r4 00809 s25 = s2 - s5 00810 s34 = s3 - s4 00811 zout ( 1, j, nout1 ) = r1 + r25 + r34 00812 r = cos2 * r25 + cos4 * r34 + r1 00813 s = sin2 * s25 + sin4 * s34 00814 zout ( 1, j, nout2 ) = r - s 00815 zout ( 1, j, nout5 ) = r + s 00816 r = cos4 * r25 + cos2 * r34 + r1 00817 s = sin4 * s25 - sin2 * s34 00818 zout ( 1, j, nout3 ) = r - s 00819 zout ( 1, j, nout4 ) = r + s 00820 r25 = r2 - r5 00821 r34 = r3 - r4 00822 s25 = s2 + s5 00823 s34 = s3 + s4 00824 zout ( 2, j, nout1 ) = s1 + s25 + s34 00825 r = cos2 * s25 + cos4 * s34 + s1 00826 s = sin2 * r25 + sin4 * r34 00827 zout ( 2, j, nout2 ) = r + s 00828 zout ( 2, j, nout5 ) = r - s 00829 r = cos4 * s25 + cos2 * s34 + s1 00830 s = sin4 * r25 - sin2 * r34 00831 zout ( 2, j, nout3 ) = r + s 00832 zout ( 2, j, nout4 ) = r - s 00833 END DO 00834 END DO 00835 ELSE 00836 nin1 = ia - after 00837 nout1 = ia - atn 00838 DO ib = 1, before 00839 nin1 = nin1 + after 00840 nin2 = nin1 + atb 00841 nin3 = nin2 + atb 00842 nin4 = nin3 + atb 00843 nin5 = nin4 + atb 00844 nout1 = nout1 + atn 00845 nout2 = nout1 + after 00846 nout3 = nout2 + after 00847 nout4 = nout3 + after 00848 nout5 = nout4 + after 00849 DO j = 1, nfft 00850 r1 = zin ( 1, j, nin1 ) 00851 s1 = zin ( 2, j, nin1 ) 00852 r = zin ( 1, j, nin2 ) 00853 s = zin ( 2, j, nin2 ) 00854 r2 = ( r + s ) * rt2i 00855 s2 = ( -r + s ) * rt2i 00856 r3 = zin ( 2, j, nin3 ) 00857 s3 = -zin ( 1, j, nin3 ) 00858 r = zin ( 1, j, nin4 ) 00859 s = zin ( 2, j, nin4 ) 00860 r4 = ( s - r ) * rt2i 00861 s4 = - ( r + s ) * rt2i 00862 r5 = -zin ( 1, j, nin5 ) 00863 s5 = -zin ( 2, j, nin5 ) 00864 r25 = r2 + r5 00865 r34 = r3 + r4 00866 s25 = s2 - s5 00867 s34 = s3 - s4 00868 zout ( 1, j, nout1 ) = r1 + r25 + r34 00869 r = cos2 * r25 + cos4 * r34 + r1 00870 s = sin2 * s25 + sin4 * s34 00871 zout ( 1, j, nout2 ) = r - s 00872 zout ( 1, j, nout5 ) = r + s 00873 r = cos4 * r25 + cos2 * r34 + r1 00874 s = sin4 * s25 - sin2 * s34 00875 zout ( 1, j, nout3 ) = r - s 00876 zout ( 1, j, nout4 ) = r + s 00877 r25 = r2 - r5 00878 r34 = r3 - r4 00879 s25 = s2 + s5 00880 s34 = s3 + s4 00881 zout ( 2, j, nout1) = s1 + s25 + s34 00882 r = cos2 * s25 + cos4 * s34 + s1 00883 s = sin2 * r25 + sin4 * r34 00884 zout ( 2, j, nout2 ) = r + s 00885 zout ( 2, j, nout5 ) = r - s 00886 r = cos4 * s25 + cos2 * s34 + s1 00887 s = sin4 * r25 - sin2 * r34 00888 zout ( 2, j, nout3 ) = r + s 00889 zout ( 2, j, nout4 ) = r - s 00890 END DO 00891 END DO 00892 END IF 00893 ELSE 00894 ias = ia - 1 00895 itt = ias * before 00896 itrig = itt + 1 00897 cr2 = trig ( 1, itrig ) 00898 ci2 = trig ( 2, itrig ) 00899 itrig = itrig + itt 00900 cr3 = trig ( 1, itrig ) 00901 ci3 = trig ( 2, itrig ) 00902 itrig = itrig + itt 00903 cr4 = trig ( 1, itrig ) 00904 ci4 = trig ( 2, itrig ) 00905 itrig = itrig + itt 00906 cr5 = trig ( 1, itrig ) 00907 ci5 = trig ( 2, itrig ) 00908 nin1 = ia - after 00909 nout1 = ia - atn 00910 DO ib = 1, before 00911 nin1 = nin1 + after 00912 nin2 = nin1 + atb 00913 nin3 = nin2 + atb 00914 nin4 = nin3 + atb 00915 nin5 = nin4 + atb 00916 nout1 = nout1 + atn 00917 nout2 = nout1 + after 00918 nout3 = nout2 + after 00919 nout4 = nout3 + after 00920 nout5 = nout4 + after 00921 DO j = 1, nfft 00922 r1 = zin ( 1, j, nin1 ) 00923 s1 = zin ( 2, j, nin1 ) 00924 r = zin ( 1, j, nin2 ) 00925 s = zin ( 2, j, nin2 ) 00926 r2 = r * cr2 - s * ci2 00927 s2 = r * ci2 + s * cr2 00928 r = zin ( 1, j, nin3 ) 00929 s = zin ( 2, j, nin3 ) 00930 r3 = r * cr3 - s * ci3 00931 s3 = r * ci3 + s * cr3 00932 r = zin ( 1, j, nin4 ) 00933 s = zin ( 2, j, nin4 ) 00934 r4 = r * cr4 - s * ci4 00935 s4 = r * ci4 + s * cr4 00936 r = zin ( 1, j, nin5 ) 00937 s = zin ( 2, j, nin5 ) 00938 r5 = r * cr5 - s * ci5 00939 s5 = r * ci5 + s * cr5 00940 r25 = r2 + r5 00941 r34 = r3 + r4 00942 s25 = s2 - s5 00943 s34 = s3 - s4 00944 zout ( 1, j, nout1 ) = r1 + r25 + r34 00945 r = cos2 * r25 + cos4 * r34 + r1 00946 s = sin2 * s25 + sin4 * s34 00947 zout ( 1, j, nout2 ) = r - s 00948 zout ( 1, j, nout5 ) = r + s 00949 r = cos4 * r25 + cos2 * r34 + r1 00950 s = sin4 * s25 - sin2 * s34 00951 zout ( 1, j, nout3 ) = r - s 00952 zout ( 1, j, nout4 ) = r + s 00953 r25 = r2 - r5 00954 r34 = r3 - r4 00955 s25 = s2 + s5 00956 s34 = s3 + s4 00957 zout ( 2, j, nout1 ) = s1 + s25 + s34 00958 r = cos2 * s25 + cos4 * s34 + s1 00959 s = sin2 * r25 + sin4 * r34 00960 zout ( 2, j, nout2 ) = r + s 00961 zout ( 2, j, nout5 ) = r - s 00962 r = cos4 * s25 + cos2 * s34 + s1 00963 s = sin4 * r25 - sin2 * r34 00964 zout ( 2, j, nout3 ) = r + s 00965 zout ( 2, j, nout4 ) = r - s 00966 END DO 00967 END DO 00968 END IF 00969 END DO 00970 ELSE IF ( now == 6 ) THEN 00971 bbs = isign * bb 00972 ia = 1 00973 nin1 = ia - after 00974 nout1 = ia - atn 00975 DO ib = 1, before 00976 nin1 = nin1 + after 00977 nin2 = nin1 + atb 00978 nin3 = nin2 + atb 00979 nin4 = nin3 + atb 00980 nin5 = nin4 + atb 00981 nin6 = nin5 + atb 00982 nout1 = nout1 + atn 00983 nout2 = nout1 + after 00984 nout3 = nout2 + after 00985 nout4 = nout3 + after 00986 nout5 = nout4 + after 00987 nout6 = nout5 + after 00988 DO j = 1, nfft 00989 r2 = zin ( 1, j, nin3 ) 00990 s2 = zin ( 2, j, nin3 ) 00991 r3 = zin ( 1, j, nin5 ) 00992 s3 = zin ( 2, j, nin5 ) 00993 r = r2 + r3 00994 s = s2 + s3 00995 r1 = zin ( 1, j, nin1 ) 00996 s1 = zin ( 2, j, nin1 ) 00997 ur1 = r + r1 00998 ui1 = s + s1 00999 r1 = r1 - 0.5_dp * r 01000 s1 = s1 - 0.5_dp * s 01001 r = r2 - r3 01002 s = s2 - s3 01003 ur2 = r1 - s * bbs 01004 ui2 = s1 + r * bbs 01005 ur3 = r1 + s * bbs 01006 ui3 = s1 - r * bbs 01007 01008 r2 = zin ( 1, j, nin6 ) 01009 s2 = zin ( 2, j, nin6 ) 01010 r3 = zin ( 1, j, nin2 ) 01011 s3 = zin ( 2, j, nin2 ) 01012 r = r2 + r3 01013 s = s2 + s3 01014 r1 = zin ( 1, j, nin4 ) 01015 s1 = zin ( 2, j, nin4 ) 01016 vr1 = r + r1 01017 vi1 = s + s1 01018 r1 = r1 - 0.5_dp * r 01019 s1 = s1 - 0.5_dp * s 01020 r = r2 - r3 01021 s = s2 - s3 01022 vr2 = r1 - s * bbs 01023 vi2 = s1 + r * bbs 01024 vr3 = r1 + s * bbs 01025 vi3 = s1 - r * bbs 01026 01027 zout ( 1, j, nout1 ) = ur1 + vr1 01028 zout ( 2, j, nout1 ) = ui1 + vi1 01029 zout ( 1, j, nout5 ) = ur2 + vr2 01030 zout ( 2, j, nout5 ) = ui2 + vi2 01031 zout ( 1, j, nout3 ) = ur3 + vr3 01032 zout ( 2, j, nout3 ) = ui3 + vi3 01033 zout ( 1, j, nout4 ) = ur1 - vr1 01034 zout ( 2, j, nout4 ) = ui1 - vi1 01035 zout ( 1, j, nout2 ) = ur2 - vr2 01036 zout ( 2, j, nout2 ) = ui2 - vi2 01037 zout ( 1, j, nout6 ) = ur3 - vr3 01038 zout ( 2, j, nout6 ) = ui3 - vi3 01039 END DO 01040 END DO 01041 ELSE 01042 STOP 'Error fftstp' 01043 END IF 01044 01045 !-----------------------------------------------------------------------------! 01046 01047 END SUBROUTINE fftstp 01048 01049 !-----------------------------------------------------------------------------!
1.7.3