Changeset 3725
- Timestamp:
- Feb 7, 2019 10:11:02 AM (5 years ago)
- Location:
- palm/trunk
- Files:
-
- 8 edited
-
SCRIPTS/palmbuild (modified) (3 diffs)
-
SCRIPTS/palmrun (modified) (3 diffs)
-
SOURCE/buoyancy.f90 (modified) (3 diffs)
-
SOURCE/data_log.f90 (modified) (7 diffs)
-
SOURCE/gust_mod.f90 (modified) (2 diffs)
-
SOURCE/poismg_mod.f90 (modified) (3 diffs)
-
SOURCE/temperton_fft_mod.f90 (modified) (23 diffs)
-
SOURCE/wind_turbine_model_mod.f90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SCRIPTS/palmbuild
r3549 r3725 27 27 # ----------------- 28 28 # $Id$ 29 # ssh-calls for compilations on remote systems modified to avoid output 30 # of login messages on specific systems changed again (reverted as before r3549) 31 # 32 # 3549 2018-11-21 15:44:44Z raasch 29 33 # ssh-calls for compilations on remote systems modified to avoid output 30 34 # of login messages on specific systems … … 571 575 fi 572 576 make_call_string="make -f Makefile_utilities $make_options F90=$compiler_name F90_SER=$compiler_name_ser COPT=\"$cpp_options\" F90FLAGS=\"$compiler_options\" LDFLAGS=\"$linker_options\" " 573 echo "$login_init_cmd $module_commands cd ${make_depository}; echo $make_call_string > LAST_MAKE_CALL; chmod u+x LAST_MAKE_CALL; $make_call_string; [[ \$? != 0 ]] && echo MAKE_ERROR" | ssh -q $ssh_key ${remote_username}@${remote_ip} 2>/dev/null | tee ${configuration_identifier}_last_make_protocol574 ###ssh -q $ssh_key ${remote_username}@${remote_ip} "$login_init_cmd $module_commands cd ${make_depository}; echo $make_call_string > LAST_MAKE_CALL; chmod u+x LAST_MAKE_CALL; $make_call_string; [[ \$? != 0 ]] && echo MAKE_ERROR" 2>&1 | tee ${configuration_identifier}_last_make_protocol577 ### echo "$login_init_cmd $module_commands cd ${make_depository}; echo $make_call_string > LAST_MAKE_CALL; chmod u+x LAST_MAKE_CALL; $make_call_string; [[ \$? != 0 ]] && echo MAKE_ERROR" | ssh -q $ssh_key ${remote_username}@${remote_ip} 2>/dev/null | tee ${configuration_identifier}_last_make_protocol 578 ssh -q $ssh_key ${remote_username}@${remote_ip} "$login_init_cmd $module_commands cd ${make_depository}; echo $make_call_string > LAST_MAKE_CALL; chmod u+x LAST_MAKE_CALL; $make_call_string; [[ \$? != 0 ]] && echo MAKE_ERROR" 2>&1 | tee ${configuration_identifier}_last_make_protocol 575 579 576 580 if [[ $(grep -c MAKE_ERROR ${configuration_identifier}_last_make_protocol) != 0 ]] … … 611 615 fi 612 616 make_call_string="make $make_options PROG=$program_name F90=$compiler_name COPT=\"$cpp_options\" F90FLAGS=\"$compiler_options\" LDFLAGS=\"$linker_options\" " 613 echo "$login_init_cmd $module_commands cd ${make_depository}; echo $make_call_string > LAST_MAKE_CALL; chmod u+x LAST_MAKE_CALL; $make_call_string; [[ \$? != 0 ]] && echo MAKE_ERROR" | ssh -q $ssh_key ${remote_username}@${remote_ip} 2>&1 | tee ${configuration_identifier}_last_make_protocol614 ###ssh -q $ssh_key ${remote_username}@${remote_ip} "$login_init_cmd $module_commands cd ${make_depository}; echo $make_call_string > LAST_MAKE_CALL; chmod u+x LAST_MAKE_CALL; $make_call_string; [[ \$? != 0 ]] && echo MAKE_ERROR" 2>&1 | tee ${configuration_identifier}_last_make_protocol617 ### echo "$login_init_cmd $module_commands cd ${make_depository}; echo $make_call_string > LAST_MAKE_CALL; chmod u+x LAST_MAKE_CALL; $make_call_string; [[ \$? != 0 ]] && echo MAKE_ERROR" | ssh -q $ssh_key ${remote_username}@${remote_ip} 2>&1 | tee ${configuration_identifier}_last_make_protocol 618 ssh -q $ssh_key ${remote_username}@${remote_ip} "$login_init_cmd $module_commands cd ${make_depository}; echo $make_call_string > LAST_MAKE_CALL; chmod u+x LAST_MAKE_CALL; $make_call_string; [[ \$? != 0 ]] && echo MAKE_ERROR" 2>&1 | tee ${configuration_identifier}_last_make_protocol 615 619 616 620 if [[ $(grep -c MAKE_ERROR ${configuration_identifier}_last_make_protocol) != 0 ]] -
palm/trunk/SCRIPTS/palmrun
r3682 r3725 27 27 # ----------------- 28 28 # $Id$ 29 # error messages for failed restarts extended 30 # 31 # 3682 2019-01-18 17:01:54Z knoop 29 32 # ssh-call for submitting batch jobs on remote systems modified again to avoid output 30 33 # of login messages on specific systems … … 2799 2802 2800 2803 # CHECK, IF RESTART JOB HAS BEEN STARTED 2801 if [[ $( grep -c "+++ palmrun killed" palmrun_restart.log ) != 0 ]]2804 if [[ $( grep -c "+++ palmrun crashed" palmrun_restart.log ) != 0 ]] 2802 2805 then 2803 2806 printf "\n$dashes\n +++ creating restart run failed \n" … … 2805 2808 rm palmrun_restart.log 2806 2809 exit 2810 elif [[ $( grep -c "*** palmrun finished" palmrun_restart.log ) != 1 ]] 2811 then 2812 printf "\n$dashes\n +++ creating restart run failed, probably due to network problems\n" 2813 locat=create_restart 2814 rm palmrun_restart.log 2815 exit 2807 2816 else 2808 2817 printf "\n$dashes\n *** restart run initiated \n" -
palm/trunk/SOURCE/buoyancy.f90
r3655 r3725 25 25 ! ----------------- 26 26 ! $Id$ 27 ! unused variables removed 28 ! 29 ! 3655 2019-01-07 16:51:22Z knoop 27 30 ! nopointer option removed 28 31 ! … … 160 163 161 164 USE indices, & 162 ONLY: nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nzb, & 163 nzt 165 ONLY: nxl, nxlu, nxr, nyn, nys, nzb, nzt 164 166 165 167 … … 251 253 252 254 USE indices, & 253 ONLY: n xlg, nxrg, nyng, nysg, nzb, nzt255 ONLY: nzb, nzt 254 256 255 257 IMPLICIT NONE -
palm/trunk/SOURCE/data_log.f90
r3655 r3725 25 25 ! ----------------- 26 26 ! $Id$ 27 ! preprocessor directives removed to avoid compiler warnings 28 ! 29 ! 3655 2019-01-07 16:51:22Z knoop 27 30 ! Corrected "Former revisions" section 28 31 ! … … 61 64 SUBROUTINE data_log( array, i1, i2, j1, j2, k1, k2 ) 62 65 63 #if defined( __logging )64 65 66 USE control_parameters, & 66 67 ONLY: log_message, simulated_time … … 98 99 WRITE ( 20 ) array 99 100 100 #endif101 101 END SUBROUTINE data_log 102 102 … … 110 110 111 111 SUBROUTINE data_log_2d( array, i1, i2, j1, j2) 112 113 #if defined( __logging )114 112 115 113 USE control_parameters, & … … 146 144 WRITE ( 20 ) array 147 145 148 #endif149 146 END SUBROUTINE data_log_2d 150 147 … … 158 155 159 156 SUBROUTINE data_log_2d_int( array, i1, i2, j1, j2) 160 161 #if defined( __logging )162 157 163 158 USE control_parameters, & … … 194 189 WRITE ( 20 ) array 195 190 196 #endif197 191 END SUBROUTINE data_log_2d_int -
palm/trunk/SOURCE/gust_mod.f90
r3685 r3725 25 25 ! ----------------- 26 26 ! $Id$ 27 ! dummy statement modified to avoid compiler warnings about unused variables 28 ! 29 ! 3685 2019-01-21 01:02:11Z knoop 27 30 ! Some interface calls moved to module_interface + cleanup 28 31 ! … … 445 448 ! 446 449 !-- Next line is just to avoid compiler warnings about unused variables. You may remove it. 447 IF ( found .AND. two_d ) idum = av + LEN( variable ) + LEN( grid ) + local_pf(nxl,nys,nzb_do) + fill_value 450 IF ( found .AND. two_d ) THEN 451 idum = av + LEN( variable ) + LEN( grid // mode ) + local_pf(nxl,nys,nzb_do) + fill_value 452 ENDIF 448 453 449 454 END SUBROUTINE gust_data_output_2d -
palm/trunk/SOURCE/poismg_mod.f90
r3655 r3725 25 25 ! ----------------- 26 26 ! $Id$ 27 ! unused subroutine removed 28 ! 29 ! 3655 2019-01-07 16:51:22Z knoop 27 30 ! unnecessary check eliminated 28 31 ! … … 143 146 INTERFACE sort_k_to_even_odd_blocks 144 147 MODULE PROCEDURE sort_k_to_even_odd_blocks 145 MODULE PROCEDURE sort_k_to_even_odd_blocks_int148 ! MODULE PROCEDURE sort_k_to_even_odd_blocks_int 146 149 MODULE PROCEDURE sort_k_to_even_odd_blocks_1d 147 150 END INTERFACE sort_k_to_even_odd_blocks … … 1152 1155 !> Version for 2D-INTEGER arrays 1153 1156 !------------------------------------------------------------------------------! 1154 SUBROUTINE sort_k_to_even_odd_blocks_int( i_mg , glevel )1155 1156 1157 USE indices, &1158 ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg1159 1160 IMPLICIT NONE1161 1162 INTEGER(iwp), INTENT(IN) :: glevel !< grid level1163 1164 INTEGER(iwp), DIMENSION(nzb:nzt_mg(glevel)+1, &1165 nys_mg(glevel)-1:nyn_mg(glevel)+1, &1166 nxl_mg(glevel)-1:nxr_mg(glevel)+1) :: &1167 i_mg !< array to be sorted 1157 ! SUBROUTINE sort_k_to_even_odd_blocks_int( i_mg , glevel ) 1158 ! 1159 ! 1160 ! USE indices, & 1161 ! ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg 1162 ! 1163 ! IMPLICIT NONE 1164 ! 1165 ! INTEGER(iwp), INTENT(IN) :: glevel !< grid level 1166 ! 1167 ! INTEGER(iwp), DIMENSION(nzb:nzt_mg(glevel)+1, & 1168 ! nys_mg(glevel)-1:nyn_mg(glevel)+1, & 1169 ! nxl_mg(glevel)-1:nxr_mg(glevel)+1) :: & 1170 ! i_mg !< array to be sorted 1168 1171 ! 1169 1172 !-- Local variables 1170 INTEGER(iwp) :: i !< index variabel along x1171 INTEGER(iwp) :: j !< index variable along y1172 INTEGER(iwp) :: k !< index variable along z1173 INTEGER(iwp) :: l !< grid level 1174 INTEGER(iwp) :: ind !< index variable along z1175 INTEGER(iwp),DIMENSION(nzb:nzt_mg(glevel)+1) :: tmp !< temporary odd-even sorted array1176 1177 1178 CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' )1179 1180 l = glevel1181 ind_even_odd = even_odd_level(l)1182 1183 DO i = nxl_mg(l)-1, nxr_mg(l)+11184 DO j = nys_mg(l)-1, nyn_mg(l)+11185 1173 ! INTEGER(iwp) :: i !< index variabel along x 1174 ! INTEGER(iwp) :: j !< index variable along y 1175 ! INTEGER(iwp) :: k !< index variable along z 1176 ! INTEGER(iwp) :: l !< grid level 1177 ! INTEGER(iwp) :: ind !< index variable along z 1178 ! INTEGER(iwp),DIMENSION(nzb:nzt_mg(glevel)+1) :: tmp !< temporary odd-even sorted array 1179 ! 1180 ! 1181 ! CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' ) 1182 ! 1183 ! l = glevel 1184 ! ind_even_odd = even_odd_level(l) 1185 ! 1186 ! DO i = nxl_mg(l)-1, nxr_mg(l)+1 1187 ! DO j = nys_mg(l)-1, nyn_mg(l)+1 1188 ! 1186 1189 ! 1187 1190 !-- Sort the data with even k index 1188 ind = nzb-11189 DO k = nzb, nzt_mg(l), 21190 ind = ind + 11191 tmp(ind) = i_mg(k,j,i)1192 ENDDO1191 ! ind = nzb-1 1192 ! DO k = nzb, nzt_mg(l), 2 1193 ! ind = ind + 1 1194 ! tmp(ind) = i_mg(k,j,i) 1195 ! ENDDO 1193 1196 ! 1194 1197 !-- Sort the data with odd k index 1195 DO k = nzb+1, nzt_mg(l)+1, 21196 ind = ind + 11197 tmp(ind) = i_mg(k,j,i)1198 ENDDO1199 1200 i_mg(:,j,i) = tmp1201 1202 ENDDO1203 ENDDO1204 1205 CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' )1206 1207 END SUBROUTINE sort_k_to_even_odd_blocks_int1198 ! DO k = nzb+1, nzt_mg(l)+1, 2 1199 ! ind = ind + 1 1200 ! tmp(ind) = i_mg(k,j,i) 1201 ! ENDDO 1202 ! 1203 ! i_mg(:,j,i) = tmp 1204 ! 1205 ! ENDDO 1206 ! ENDDO 1207 ! 1208 ! CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' ) 1209 ! 1210 ! END SUBROUTINE sort_k_to_even_odd_blocks_int 1208 1211 1209 1212 -
palm/trunk/SOURCE/temperton_fft_mod.f90
r3241 r3725 9 9 ! ----------------- 10 10 ! $Id$ 11 ! GOTO statements replaced, file re-formatted corresponding to coding standards 12 ! 13 ! 3241 2018-09-12 15:02:00Z raasch 11 14 ! unused variables removed 12 15 ! … … 25 28 ! Module renamed 26 29 ! 27 !28 30 ! 1682 2015-10-07 23:56:08Z knoop 29 31 ! Code annotations made doxygen readable … … 70 72 ! Description: 71 73 ! ------------ 72 !> Calls fortran-versions of fft's.74 !> Calls Fortran-versions of fft's. 73 75 !> 74 76 !> Method: … … 115 117 !> dimension a(n),work(n),trigs(n),ifax(1) 116 118 !------------------------------------------------------------------------------! 117 SUBROUTINE fft991cy(a,work,trigs,ifax,inc,jump,n,lot,isign)119 SUBROUTINE fft991cy( a, work, trigs, ifax, inc, jump, n, lot, isign ) 118 120 119 121 USE kinds … … 121 123 IMPLICIT NONE 122 124 123 ! Scalar arguments124 125 INTEGER(iwp) :: inc !< 125 126 INTEGER(iwp) :: isign !< … … 128 129 INTEGER(iwp) :: n !< 129 130 130 ! Array arguments131 131 REAL(wp) :: a(*) !< 132 132 REAL(wp) :: trigs(*) !< … … 134 134 INTEGER(iwp) :: ifax(*) !< 135 135 136 ! Local scalars:137 136 INTEGER(iwp) :: i !< 138 137 INTEGER(iwp) :: ia !< … … 156 155 INTEGER(iwp) :: nx !< 157 156 158 ! Intrinsic functions 159 ! INTRINSIC MOD 160 161 162 ! Executable statements 163 164 IF (ifax(10)/=n) CALL set99(trigs,ifax,n) 157 158 IF ( ifax(10) /= n ) CALL set99( trigs, ifax, n ) 165 159 nfax = ifax(1) 166 nx = n + 1 167 IF (MOD(n,2)==1) nx = n 168 nblox = 1 + (lot-1)/nfft 169 nvex = lot - (nblox-1)*nfft 170 IF (isign==-1) GO TO 50 171 172 ! isign=+1, spectral to gridpoint transform 173 174 istart = 1 175 DO nb = 1, nblox 176 ia = istart 177 i = istart 178 !DIR$ IVDEP 179 !CDIR NODEP 180 !OCL NOVREC 181 DO j = 1, nvex 182 a(i+inc) = 0.5_wp*a(i) 183 i = i + jump 184 END DO 185 IF (MOD(n,2)==1) GO TO 10 186 i = istart + n*inc 187 DO j = 1, nvex 188 a(i) = 0.5_wp*a(i) 189 i = i + jump 190 END DO 191 10 CONTINUE 192 ia = istart + inc 193 la = 1 194 igo = + 1 195 196 DO k = 1, nfax 197 ifac = ifax(k+1) 198 ierr = -1 199 IF (igo==-1) GO TO 20 200 CALL rpassm(a(ia),a(ia+la*inc),work(1),work(ifac*la+1),trigs,inc,1, & 201 & jump,nx,nvex,n,ifac,la,ierr) 202 GO TO 30 203 20 CONTINUE 204 CALL rpassm(work(1),work(la+1),a(ia),a(ia+ifac*la*inc),trigs,1,inc,nx, & 205 & jump,nvex,n,ifac,la,ierr) 206 30 CONTINUE 207 IF (ierr/=0) GO TO 100 208 la = ifac*la 209 igo = -igo 160 nx = n + 1 161 IF ( MOD(n,2) == 1 ) nx = n 162 nblox = 1 + (lot-1) / nfft 163 nvex = lot - (nblox-1) * nfft 164 165 166 IF ( isign == 1 ) THEN 167 ! 168 !-- Backward fft: spectral to gridpoint transform 169 170 istart = 1 171 DO nb = 1, nblox 172 210 173 ia = istart 211 END DO 212 213 ! If necessary, copy results back to a 214 215 IF (MOD(nfax,2)==0) GO TO 40 216 ibase = 1 217 jbase = ia 218 DO jj = 1, nvex 219 i = ibase 220 j = jbase 221 DO ii = 1, n 222 a(j) = work(i) 223 i = i + 1 224 j = j + inc 225 END DO 226 ibase = ibase + nx 227 jbase = jbase + jump 228 END DO 229 40 CONTINUE 230 231 ! Fill in zeros at end 232 233 ix = istart + n*inc 234 !DIR$ IVDEP 235 !CDIR NODEP 236 !OCL NOVREC 237 DO j = 1, nvex 238 a(ix) = 0.0_wp 239 a(ix+inc) = 0.0_wp 240 ix = ix + jump 241 END DO 242 243 istart = istart + nvex*jump 244 nvex = nfft 245 END DO 246 RETURN 247 248 ! isign=-1, gridpoint to spectral transform 249 250 50 CONTINUE 251 istart = 1 252 DO nb = 1, nblox 253 ia = istart 254 la = n 255 igo = + 1 256 257 DO k = 1, nfax 258 ifac = ifax(nfax+2-k) 259 la = la/ifac 260 ierr = -1 261 IF (igo==-1) GO TO 60 262 CALL qpassm(a(ia),a(ia+ifac*la*inc),work(1),work(la+1),trigs,inc,1, & 263 & jump,nx,nvex,n,ifac,la,ierr) 264 GO TO 70 265 60 CONTINUE 266 CALL qpassm(work(1),work(ifac*la+1),a(ia),a(ia+la*inc),trigs,1,inc,nx, & 267 & jump,nvex,n,ifac,la,ierr) 268 70 CONTINUE 269 IF (ierr/=0) GO TO 100 270 igo = -igo 174 i = istart 175 !DIR$ IVDEP 176 !CDIR NODEP 177 !OCL NOVREC 178 DO j = 1, nvex 179 a(i+inc) = 0.5_wp * a(i) 180 i = i + jump 181 ENDDO 182 IF ( MOD(n,2) /= 1 ) THEN 183 i = istart + n * inc 184 DO j = 1, nvex 185 a(i) = 0.5_wp * a(i) 186 i = i + jump 187 ENDDO 188 ENDIF 271 189 ia = istart + inc 272 END DO 273 274 ! If necessary, copy results back to a 275 276 IF (MOD(nfax,2)==0) GO TO 80 277 ibase = 1 278 jbase = ia 279 DO jj = 1, nvex 280 i = ibase 281 j = jbase 282 DO ii = 1, n 283 a(j) = work(i) 284 i = i + 1 285 j = j + inc 286 END DO 287 ibase = ibase + nx 288 jbase = jbase + jump 289 END DO 290 80 CONTINUE 291 292 ! Shift a(0) & fill in zero imag parts 293 294 ix = istart 295 !DIR$ IVDEP 296 !CDIR NODEP 297 !OCL NOVREC 298 DO j = 1, nvex 299 a(ix) = a(ix+inc) 300 a(ix+inc) = 0.0_wp 301 ix = ix + jump 302 END DO 303 IF (MOD(n,2)==1) GO TO 90 304 iz = istart + (n+1)*inc 305 DO j = 1, nvex 306 a(iz) = 0.0_wp 307 iz = iz + jump 308 END DO 309 90 CONTINUE 310 311 istart = istart + nvex*jump 312 nvex = nfft 313 END DO 314 RETURN 315 316 ! Error messages 317 318 100 CONTINUE 319 320 SELECT CASE (ierr) 321 CASE (:-1) 322 WRITE (nout,'(A,I5,A)') ' Vector length =',nvex,', greater than nfft' 323 CASE (0) 324 WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', not catered for' 325 CASE (1:) 326 WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', only catered for if la*ifac=n' 327 END SELECT 328 329 RETURN 330 END SUBROUTINE fft991cy 190 la = 1 191 igo = + 1 192 193 DO k = 1, nfax 194 195 ifac = ifax(k+1) 196 ierr = -1 197 IF ( igo /= -1 ) THEN 198 CALL rpassm( a(ia), a(ia+la*inc), work(1), work(ifac*la+1), trigs, inc, 1, jump, & 199 nx, nvex, n, ifac, la, ierr ) 200 ELSE 201 CALL rpassm( work(1), work(la+1), a(ia), a(ia+ifac*la*inc), trigs, 1, inc, nx, & 202 jump, nvex, n, ifac, la, ierr ) 203 ENDIF 204 ! 205 !-- Following messages shouldn't appear in PALM applications 206 IF ( ierr /= 0 ) THEN 207 SELECT CASE (ierr) 208 CASE (:-1) 209 WRITE (nout,'(A,I5,A)') ' Vector length =',nvex,', greater than nfft' 210 CASE (0) 211 WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', not catered for' 212 CASE (1:) 213 WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', only catered for if la*ifac=n' 214 END SELECT 215 RETURN 216 ENDIF 217 218 la = ifac * la 219 igo = -igo 220 ia = istart 221 222 ENDDO 223 ! 224 !-- If necessary, copy results back to a 225 IF ( MOD(nfax,2) /= 0 ) THEN 226 ibase = 1 227 jbase = ia 228 DO jj = 1, nvex 229 i = ibase 230 j = jbase 231 DO ii = 1, n 232 a(j) = work(i) 233 i = i + 1 234 j = j + inc 235 ENDDO 236 ibase = ibase + nx 237 jbase = jbase + jump 238 ENDDO 239 ENDIF 240 ! 241 !-- Fill in zeros at end 242 ix = istart + n*inc 243 !DIR$ IVDEP 244 !CDIR NODEP 245 !OCL NOVREC 246 DO j = 1, nvex 247 a(ix) = 0.0_wp 248 a(ix+inc) = 0.0_wp 249 ix = ix + jump 250 ENDDO 251 istart = istart + nvex*jump 252 nvex = nfft 253 254 ENDDO 255 256 ELSEIF ( isign == -1 ) THEN 257 ! 258 !-- Forward fft: gridpoint to spectral transform 259 istart = 1 260 DO nb = 1, nblox 261 ia = istart 262 la = n 263 igo = + 1 264 265 DO k = 1, nfax 266 267 ifac = ifax(nfax+2-k) 268 la = la / ifac 269 ierr = -1 270 IF ( igo /= -1 ) THEN 271 CALL qpassm( a(ia), a(ia+ifac*la*inc), work(1), work(la+1), trigs, inc, 1, jump, & 272 nx, nvex, n, ifac, la, ierr ) 273 ELSE 274 CALL qpassm( work(1), work(ifac*la+1), a(ia), a(ia+la*inc), trigs, 1, inc, nx, & 275 jump, nvex, n, ifac, la, ierr ) 276 ENDIF 277 ! 278 !-- Following messages shouldn't appear in PALM applications 279 IF ( ierr /= 0 ) THEN 280 SELECT CASE (ierr) 281 CASE (:-1) 282 WRITE (nout,'(A,I5,A)') ' Vector length =',nvex,', greater than nfft' 283 CASE (0) 284 WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', not catered for' 285 CASE (1:) 286 WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', only catered for if la*ifac=n' 287 END SELECT 288 RETURN 289 ENDIF 290 291 igo = -igo 292 ia = istart + inc 293 294 ENDDO 295 ! 296 !-- If necessary, copy results back to a 297 IF ( MOD(nfax,2) /= 0 ) THEN 298 ibase = 1 299 jbase = ia 300 DO jj = 1, nvex 301 i = ibase 302 j = jbase 303 DO ii = 1, n 304 a(j) = work(i) 305 i = i + 1 306 j = j + inc 307 ENDDO 308 ibase = ibase + nx 309 jbase = jbase + jump 310 ENDDO 311 ENDIF 312 ! 313 !-- Shift a(0) and fill in zero imag parts 314 ix = istart 315 !DIR$ IVDEP 316 !CDIR NODEP 317 !OCL NOVREC 318 DO j = 1, nvex 319 a(ix) = a(ix+inc) 320 a(ix+inc) = 0.0_wp 321 ix = ix + jump 322 ENDDO 323 IF ( MOD(n,2) /= 1 ) THEN 324 iz = istart + (n+1) * inc 325 DO j = 1, nvex 326 a(iz) = 0.0_wp 327 iz = iz + jump 328 ENDDO 329 ENDIF 330 istart = istart + nvex * jump 331 nvex = nfft 332 333 ENDDO 334 335 ENDIF 336 337 END SUBROUTINE fft991cy 331 338 332 339 !------------------------------------------------------------------------------! … … 357 364 !> 3 - ifac only catered for if la=n/ifac 358 365 !------------------------------------------------------------------------------! 359 SUBROUTINE qpassm(a,b,c,d,trigs,inc1,inc2,inc3,inc4,lot,n,ifac,la,ierr)366 SUBROUTINE qpassm( a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr ) 360 367 361 368 USE kinds … … 363 370 IMPLICIT NONE 364 371 365 ! Scalar arguments366 372 INTEGER(iwp) :: ierr !< 367 373 INTEGER(iwp) :: ifac !< … … 374 380 INTEGER(iwp) :: n !< 375 381 376 ! Array arguments 377 ! REAL :: a(n),b(n),c(n),d(n),trigs(n) 382 ! 383 !-- Arrays are dimensioned with n 378 384 REAL(wp) :: a(*) !< 379 385 REAL(wp) :: b(*) !< … … 382 388 REAL(wp) :: trigs(*) !< 383 389 384 ! Local scalars:385 390 REAL(wp) :: a0 !< 386 391 REAL(wp) :: a1 !< … … 430 435 INTEGER(iwp) :: ia !< 431 436 INTEGER(iwp) :: ib !< 432 INTEGER(iwp) :: ibad !<433 437 INTEGER(iwp) :: ibase !< 434 438 INTEGER(iwp) :: ic !< … … 461 465 INTEGER(iwp) :: m !< 462 466 463 ! Intrinsic functions 464 ! INTRINSIC REAL, SQRT 465 466 ! Data statements 467 DATA sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/, & 468 & qrt5/0.559016994374947_wp/, sin60/0.866025403784437_wp/469 470 471 ! Executable statements472 473 m = n/ifac474 iink = la*inc1 475 jink = la*inc2476 i jump = (ifac-1)*iink477 kstop = (n-ifac)/(2*ifac)478 479 ibad = 1480 IF (lot>nfft) GO TO 180 467 468 DATA sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/, & 469 qrt5/0.559016994374947_wp/, sin60/0.866025403784437_wp/ 470 471 472 ierr = 0 473 474 IF ( lot > nfft ) THEN 475 ierr = 1 476 RETURN 477 ENDIF 478 479 m = n / ifac 480 iink = la * inc1 481 jink = la * inc2 482 ijump = (ifac-1) * iink 483 kstop = ( n-ifac ) / ( 2*ifac ) 484 481 485 ibase = 0 482 486 jbase = 0 483 487 igo = ifac - 1 484 IF (igo==7) igo = 6 485 ibad = 2 486 IF (igo<1 .OR. igo>6) GO TO 180 487 GO TO (10,40,70,100,130,160) igo 488 489 ! Coding for factor 2 490 491 10 CONTINUE 492 ia = 1 493 ib = ia + iink 494 ja = 1 495 jb = ja + (2*m-la)*inc2 496 497 IF (la==m) GO TO 30 498 499 DO l = 1, la 500 i = ibase 501 j = jbase 502 !DIR$ IVDEP 503 !CDIR NODEP 504 !OCL NOVREC 505 DO ijk = 1, lot 506 c(ja+j) = a(ia+i) + a(ib+i) 507 c(jb+j) = a(ia+i) - a(ib+i) 508 i = i + inc3 509 j = j + inc4 510 END DO 511 ibase = ibase + inc1 512 jbase = jbase + inc2 513 END DO 514 ja = ja + jink 515 jink = 2*jink 516 jb = jb - jink 517 ibase = ibase + ijump 518 ijump = 2*ijump + iink 519 IF (ja==jb) GO TO 20 520 DO k = la, kstop, la 521 kb = k + k 522 c1 = trigs(kb+1) 523 s1 = trigs(kb+2) 524 jbase = 0 525 DO l = 1, la 526 i = ibase 527 j = jbase 528 !DIR$ IVDEP 529 !CDIR NODEP 530 !OCL NOVREC 531 DO ijk = 1, lot 532 c(ja+j) = a(ia+i) + (c1*a(ib+i)+s1*b(ib+i)) 533 c(jb+j) = a(ia+i) - (c1*a(ib+i)+s1*b(ib+i)) 534 d(ja+j) = (c1*b(ib+i)-s1*a(ib+i)) + b(ia+i) 535 d(jb+j) = (c1*b(ib+i)-s1*a(ib+i)) - b(ia+i) 536 i = i + inc3 537 j = j + inc4 538 END DO 539 ibase = ibase + inc1 540 jbase = jbase + inc2 541 END DO 542 ibase = ibase + ijump 543 ja = ja + jink 544 jb = jb - jink 545 END DO 546 IF (ja>jb) GO TO 170 547 20 CONTINUE 548 jbase = 0 549 DO l = 1, la 550 i = ibase 551 j = jbase 552 !DIR$ IVDEP 553 !CDIR NODEP 554 !OCL NOVREC 555 DO ijk = 1, lot 556 c(ja+j) = a(ia+i) 557 d(ja+j) = -a(ib+i) 558 i = i + inc3 559 j = j + inc4 560 END DO 561 ibase = ibase + inc1 562 jbase = jbase + inc2 563 END DO 564 565 GO TO 170 566 30 CONTINUE 567 z = 1.0_wp/REAL(n,KIND=wp) 568 DO l = 1, la 569 i = ibase 570 j = jbase 571 !DIR$ IVDEP 572 !CDIR NODEP 573 !OCL NOVREC 574 DO ijk = 1, lot 575 c(ja+j) = z*(a(ia+i)+a(ib+i)) 576 c(jb+j) = z*(a(ia+i)-a(ib+i)) 577 i = i + inc3 578 j = j + inc4 579 END DO 580 ibase = ibase + inc1 581 jbase = jbase + inc2 582 END DO 583 GO TO 170 584 585 ! Coding for factor 3 586 587 40 CONTINUE 588 ia = 1 589 ib = ia + iink 590 ic = ib + iink 591 ja = 1 592 jb = ja + (2*m-la)*inc2 593 jc = jb 594 595 IF (la==m) GO TO 60 596 597 DO l = 1, la 598 i = ibase 599 j = jbase 600 !DIR$ IVDEP 601 !CDIR NODEP 602 !OCL NOVREC 603 DO ijk = 1, lot 604 c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i)) 605 c(jb+j) = a(ia+i) - 0.5_wp*(a(ib+i)+a(ic+i)) 606 d(jb+j) = sin60*(a(ic+i)-a(ib+i)) 607 i = i + inc3 608 j = j + inc4 609 END DO 610 ibase = ibase + inc1 611 jbase = jbase + inc2 612 END DO 613 ja = ja + jink 614 jink = 2*jink 615 jb = jb + jink 616 jc = jc - jink 617 ibase = ibase + ijump 618 ijump = 2*ijump + iink 619 IF (ja==jc) GO TO 50 620 DO k = la, kstop, la 621 kb = k + k 622 kc = kb + kb 623 c1 = trigs(kb+1) 624 s1 = trigs(kb+2) 625 c2 = trigs(kc+1) 626 s2 = trigs(kc+2) 627 jbase = 0 628 DO l = 1, la 629 i = ibase 630 j = jbase 631 !DIR$ IVDEP 632 !CDIR NODEP 633 !OCL NOVREC 634 DO ijk = 1, lot 635 a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c2*a(ic+i)+s2*b(ic+i)) 636 b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c2*b(ic+i)-s2*a(ic+i)) 637 a2 = a(ia+i) - 0.5_wp*a1 638 b2 = b(ia+i) - 0.5_wp*b1 639 a3 = sin60*((c1*a(ib+i)+s1*b(ib+i))-(c2*a(ic+i)+s2*b(ic+i))) 640 b3 = sin60*((c1*b(ib+i)-s1*a(ib+i))-(c2*b(ic+i)-s2*a(ic+i))) 641 c(ja+j) = a(ia+i) + a1 642 d(ja+j) = b(ia+i) + b1 643 c(jb+j) = a2 + b3 644 d(jb+j) = b2 - a3 645 c(jc+j) = a2 - b3 646 d(jc+j) = -(b2+a3) 647 i = i + inc3 648 j = j + inc4 649 END DO 650 ibase = ibase + inc1 651 jbase = jbase + inc2 652 END DO 653 ibase = ibase + ijump 654 ja = ja + jink 655 jb = jb + jink 656 jc = jc - jink 657 END DO 658 IF (ja>jc) GO TO 170 659 50 CONTINUE 660 jbase = 0 661 DO l = 1, la 662 i = ibase 663 j = jbase 664 !DIR$ IVDEP 665 !CDIR NODEP 666 !OCL NOVREC 667 DO ijk = 1, lot 668 c(ja+j) = a(ia+i) + 0.5_wp*(a(ib+i)-a(ic+i)) 669 d(ja+j) = -sin60*(a(ib+i)+a(ic+i)) 670 c(jb+j) = a(ia+i) - (a(ib+i)-a(ic+i)) 671 i = i + inc3 672 j = j + inc4 673 END DO 674 ibase = ibase + inc1 675 jbase = jbase + inc2 676 END DO 677 678 GO TO 170 679 60 CONTINUE 680 z = 1.0_wp/REAL(n,KIND=wp) 681 zsin60 = z*sin60 682 DO l = 1, la 683 i = ibase 684 j = jbase 685 !DIR$ IVDEP 686 !CDIR NODEP 687 !OCL NOVREC 688 DO ijk = 1, lot 689 c(ja+j) = z*(a(ia+i)+(a(ib+i)+a(ic+i))) 690 c(jb+j) = z*(a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i))) 691 d(jb+j) = zsin60*(a(ic+i)-a(ib+i)) 692 i = i + inc3 693 j = j + inc4 694 END DO 695 ibase = ibase + inc1 696 jbase = jbase + inc2 697 END DO 698 GO TO 170 699 700 ! Coding for factor 4 701 702 70 CONTINUE 703 ia = 1 704 ib = ia + iink 705 ic = ib + iink 706 id = ic + iink 707 ja = 1 708 jb = ja + (2*m-la)*inc2 709 jc = jb + 2*m*inc2 710 jd = jb 711 712 IF (la==m) GO TO 90 713 714 DO l = 1, la 715 i = ibase 716 j = jbase 717 !DIR$ IVDEP 718 !CDIR NODEP 719 !OCL NOVREC 720 DO ijk = 1, lot 721 c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i)) 722 c(jc+j) = (a(ia+i)+a(ic+i)) - (a(ib+i)+a(id+i)) 723 c(jb+j) = a(ia+i) - a(ic+i) 724 d(jb+j) = a(id+i) - a(ib+i) 725 i = i + inc3 726 j = j + inc4 727 END DO 728 ibase = ibase + inc1 729 jbase = jbase + inc2 730 END DO 731 ja = ja + jink 732 jink = 2*jink 733 jb = jb + jink 734 jc = jc - jink 735 jd = jd - jink 736 ibase = ibase + ijump 737 ijump = 2*ijump + iink 738 IF (jb==jc) GO TO 80 739 DO k = la, kstop, la 740 kb = k + k 741 kc = kb + kb 742 kd = kc + kb 743 c1 = trigs(kb+1) 744 s1 = trigs(kb+2) 745 c2 = trigs(kc+1) 746 s2 = trigs(kc+2) 747 c3 = trigs(kd+1) 748 s3 = trigs(kd+2) 749 jbase = 0 750 DO l = 1, la 751 i = ibase 752 j = jbase 753 !DIR$ IVDEP 754 !CDIR NODEP 755 !OCL NOVREC 756 DO ijk = 1, lot 757 a0 = a(ia+i) + (c2*a(ic+i)+s2*b(ic+i)) 758 a2 = a(ia+i) - (c2*a(ic+i)+s2*b(ic+i)) 759 a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c3*a(id+i)+s3*b(id+i)) 760 a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c3*a(id+i)+s3*b(id+i)) 761 b0 = b(ia+i) + (c2*b(ic+i)-s2*a(ic+i)) 762 b2 = b(ia+i) - (c2*b(ic+i)-s2*a(ic+i)) 763 b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c3*b(id+i)-s3*a(id+i)) 764 b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c3*b(id+i)-s3*a(id+i)) 765 c(ja+j) = a0 + a1 766 c(jc+j) = a0 - a1 767 d(ja+j) = b0 + b1 768 d(jc+j) = b1 - b0 769 c(jb+j) = a2 + b3 770 c(jd+j) = a2 - b3 771 d(jb+j) = b2 - a3 772 d(jd+j) = -(b2+a3) 773 i = i + inc3 774 j = j + inc4 775 END DO 776 ibase = ibase + inc1 777 jbase = jbase + inc2 778 END DO 779 ibase = ibase + ijump 780 ja = ja + jink 781 jb = jb + jink 782 jc = jc - jink 783 jd = jd - jink 784 END DO 785 IF (jb>jc) GO TO 170 786 80 CONTINUE 787 sin45 = SQRT(0.5_wp) 788 jbase = 0 789 DO l = 1, la 790 i = ibase 791 j = jbase 792 !DIR$ IVDEP 793 !CDIR NODEP 794 !OCL NOVREC 795 DO ijk = 1, lot 796 c(ja+j) = a(ia+i) + sin45*(a(ib+i)-a(id+i)) 797 c(jb+j) = a(ia+i) - sin45*(a(ib+i)-a(id+i)) 798 d(ja+j) = -a(ic+i) - sin45*(a(ib+i)+a(id+i)) 799 d(jb+j) = a(ic+i) - sin45*(a(ib+i)+a(id+i)) 800 i = i + inc3 801 j = j + inc4 802 END DO 803 ibase = ibase + inc1 804 jbase = jbase + inc2 805 END DO 806 807 GO TO 170 808 90 CONTINUE 809 z = 1.0_wp/REAL(n,KIND=wp) 810 DO l = 1, la 811 i = ibase 812 j = jbase 813 !DIR$ IVDEP 814 !CDIR NODEP 815 !OCL NOVREC 816 DO ijk = 1, lot 817 c(ja+j) = z*((a(ia+i)+a(ic+i))+(a(ib+i)+a(id+i))) 818 c(jc+j) = z*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) 819 c(jb+j) = z*(a(ia+i)-a(ic+i)) 820 d(jb+j) = z*(a(id+i)-a(ib+i)) 821 i = i + inc3 822 j = j + inc4 823 END DO 824 ibase = ibase + inc1 825 jbase = jbase + inc2 826 END DO 827 GO TO 170 828 829 ! Coding for factor 5 830 831 100 CONTINUE 832 ia = 1 833 ib = ia + iink 834 ic = ib + iink 835 id = ic + iink 836 ie = id + iink 837 ja = 1 838 jb = ja + (2*m-la)*inc2 839 jc = jb + 2*m*inc2 840 jd = jc 841 je = jb 842 843 IF (la==m) GO TO 120 844 845 DO l = 1, la 846 i = ibase 847 j = jbase 848 !DIR$ IVDEP 849 !CDIR NODEP 850 !OCL NOVREC 851 DO ijk = 1, lot 852 a1 = a(ib+i) + a(ie+i) 853 a3 = a(ib+i) - a(ie+i) 854 a2 = a(ic+i) + a(id+i) 855 a4 = a(ic+i) - a(id+i) 856 a5 = a(ia+i) - 0.25_wp*(a1+a2) 857 a6 = qrt5*(a1-a2) 858 c(ja+j) = a(ia+i) + (a1+a2) 859 c(jb+j) = a5 + a6 860 c(jc+j) = a5 - a6 861 d(jb+j) = -sin72*a3 - sin36*a4 862 d(jc+j) = -sin36*a3 + sin72*a4 863 i = i + inc3 864 j = j + inc4 865 END DO 866 ibase = ibase + inc1 867 jbase = jbase + inc2 868 END DO 869 ja = ja + jink 870 jink = 2*jink 871 jb = jb + jink 872 jc = jc + jink 873 jd = jd - jink 874 je = je - jink 875 ibase = ibase + ijump 876 ijump = 2*ijump + iink 877 IF (jb==jd) GO TO 110 878 DO k = la, kstop, la 879 kb = k + k 880 kc = kb + kb 881 kd = kc + kb 882 ke = kd + kb 883 c1 = trigs(kb+1) 884 s1 = trigs(kb+2) 885 c2 = trigs(kc+1) 886 s2 = trigs(kc+2) 887 c3 = trigs(kd+1) 888 s3 = trigs(kd+2) 889 c4 = trigs(ke+1) 890 s4 = trigs(ke+2) 891 jbase = 0 892 DO l = 1, la 893 i = ibase 894 j = jbase 895 !DIR$ IVDEP 896 !CDIR NODEP 897 !OCL NOVREC 898 DO ijk = 1, lot 899 a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c4*a(ie+i)+s4*b(ie+i)) 900 a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c4*a(ie+i)+s4*b(ie+i)) 901 a2 = (c2*a(ic+i)+s2*b(ic+i)) + (c3*a(id+i)+s3*b(id+i)) 902 a4 = (c2*a(ic+i)+s2*b(ic+i)) - (c3*a(id+i)+s3*b(id+i)) 903 b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c4*b(ie+i)-s4*a(ie+i)) 904 b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c4*b(ie+i)-s4*a(ie+i)) 905 b2 = (c2*b(ic+i)-s2*a(ic+i)) + (c3*b(id+i)-s3*a(id+i)) 906 b4 = (c2*b(ic+i)-s2*a(ic+i)) - (c3*b(id+i)-s3*a(id+i)) 907 a5 = a(ia+i) - 0.25_wp*(a1+a2) 908 a6 = qrt5*(a1-a2) 909 b5 = b(ia+i) - 0.25_wp*(b1+b2) 910 b6 = qrt5*(b1-b2) 911 a10 = a5 + a6 912 a20 = a5 - a6 913 b10 = b5 + b6 914 b20 = b5 - b6 915 a11 = sin72*b3 + sin36*b4 916 a21 = sin36*b3 - sin72*b4 917 b11 = sin72*a3 + sin36*a4 918 b21 = sin36*a3 - sin72*a4 919 c(ja+j) = a(ia+i) + (a1+a2) 920 c(jb+j) = a10 + a11 921 c(je+j) = a10 - a11 922 c(jc+j) = a20 + a21 923 c(jd+j) = a20 - a21 924 d(ja+j) = b(ia+i) + (b1+b2) 925 d(jb+j) = b10 - b11 926 d(je+j) = -(b10+b11) 927 d(jc+j) = b20 - b21 928 d(jd+j) = -(b20+b21) 929 i = i + inc3 930 j = j + inc4 931 END DO 932 ibase = ibase + inc1 933 jbase = jbase + inc2 934 END DO 935 ibase = ibase + ijump 936 ja = ja + jink 937 jb = jb + jink 938 jc = jc + jink 939 jd = jd - jink 940 je = je - jink 941 END DO 942 IF (jb>jd) GO TO 170 943 110 CONTINUE 944 jbase = 0 945 DO l = 1, la 946 i = ibase 947 j = jbase 948 !DIR$ IVDEP 949 !CDIR NODEP 950 !OCL NOVREC 951 DO ijk = 1, lot 952 a1 = a(ib+i) + a(ie+i) 953 a3 = a(ib+i) - a(ie+i) 954 a2 = a(ic+i) + a(id+i) 955 a4 = a(ic+i) - a(id+i) 956 a5 = a(ia+i) + 0.25_wp*(a3-a4) 957 a6 = qrt5*(a3+a4) 958 c(ja+j) = a5 + a6 959 c(jb+j) = a5 - a6 960 c(jc+j) = a(ia+i) - (a3-a4) 961 d(ja+j) = -sin36*a1 - sin72*a2 962 d(jb+j) = -sin72*a1 + sin36*a2 963 i = i + inc3 964 j = j + inc4 965 END DO 966 ibase = ibase + inc1 967 jbase = jbase + inc2 968 END DO 969 970 GO TO 170 971 120 CONTINUE 972 z = 1.0_wp/REAL(n,KIND=wp) 973 zqrt5 = z*qrt5 974 zsin36 = z*sin36 975 zsin72 = z*sin72 976 DO l = 1, la 977 i = ibase 978 j = jbase 979 !DIR$ IVDEP 980 !CDIR NODEP 981 !OCL NOVREC 982 DO ijk = 1, lot 983 a1 = a(ib+i) + a(ie+i) 984 a3 = a(ib+i) - a(ie+i) 985 a2 = a(ic+i) + a(id+i) 986 a4 = a(ic+i) - a(id+i) 987 a5 = z*(a(ia+i)-0.25_wp*(a1+a2)) 988 a6 = zqrt5*(a1-a2) 989 c(ja+j) = z*(a(ia+i)+(a1+a2)) 990 c(jb+j) = a5 + a6 991 c(jc+j) = a5 - a6 992 d(jb+j) = -zsin72*a3 - zsin36*a4 993 d(jc+j) = -zsin36*a3 + zsin72*a4 994 i = i + inc3 995 j = j + inc4 996 END DO 997 ibase = ibase + inc1 998 jbase = jbase + inc2 999 END DO 1000 GO TO 170 1001 1002 ! Coding for factor 6 1003 1004 130 CONTINUE 1005 ia = 1 1006 ib = ia + iink 1007 ic = ib + iink 1008 id = ic + iink 1009 ie = id + iink 1010 if = ie + iink 1011 ja = 1 1012 jb = ja + (2*m-la)*inc2 1013 jc = jb + 2*m*inc2 1014 jd = jc + 2*m*inc2 1015 je = jc 1016 jf = jb 1017 1018 IF (la==m) GO TO 150 1019 1020 DO l = 1, la 1021 i = ibase 1022 j = jbase 1023 !DIR$ IVDEP 1024 !CDIR NODEP 1025 !OCL NOVREC 1026 DO ijk = 1, lot 1027 a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i)) 1028 c(ja+j) = (a(ia+i)+a(id+i)) + a11 1029 c(jc+j) = (a(ia+i)+a(id+i)-0.5_wp*a11) 1030 d(jc+j) = sin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i))) 1031 a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i)) 1032 c(jb+j) = (a(ia+i)-a(id+i)) - 0.5_wp*a11 1033 d(jb+j) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i))) 1034 c(jd+j) = (a(ia+i)-a(id+i)) + a11 1035 i = i + inc3 1036 j = j + inc4 1037 END DO 1038 ibase = ibase + inc1 1039 jbase = jbase + inc2 1040 END DO 1041 ja = ja + jink 1042 jink = 2*jink 1043 jb = jb + jink 1044 jc = jc + jink 1045 jd = jd - jink 1046 je = je - jink 1047 jf = jf - jink 1048 ibase = ibase + ijump 1049 ijump = 2*ijump + iink 1050 IF (jc==jd) GO TO 140 1051 DO k = la, kstop, la 1052 kb = k + k 1053 kc = kb + kb 1054 kd = kc + kb 1055 ke = kd + kb 1056 kf = ke + kb 1057 c1 = trigs(kb+1) 1058 s1 = trigs(kb+2) 1059 c2 = trigs(kc+1) 1060 s2 = trigs(kc+2) 1061 c3 = trigs(kd+1) 1062 s3 = trigs(kd+2) 1063 c4 = trigs(ke+1) 1064 s4 = trigs(ke+2) 1065 c5 = trigs(kf+1) 1066 s5 = trigs(kf+2) 1067 jbase = 0 1068 DO l = 1, la 1069 i = ibase 1070 j = jbase 1071 !DIR$ IVDEP 1072 !CDIR NODEP 1073 !OCL NOVREC 1074 DO ijk = 1, lot 1075 a1 = c1*a(ib+i) + s1*b(ib+i) 1076 b1 = c1*b(ib+i) - s1*a(ib+i) 1077 a2 = c2*a(ic+i) + s2*b(ic+i) 1078 b2 = c2*b(ic+i) - s2*a(ic+i) 1079 a3 = c3*a(id+i) + s3*b(id+i) 1080 b3 = c3*b(id+i) - s3*a(id+i) 1081 a4 = c4*a(ie+i) + s4*b(ie+i) 1082 b4 = c4*b(ie+i) - s4*a(ie+i) 1083 a5 = c5*a(if+i) + s5*b(if+i) 1084 b5 = c5*b(if+i) - s5*a(if+i) 1085 a11 = (a2+a5) + (a1+a4) 1086 a20 = (a(ia+i)+a3) - 0.5_wp*a11 1087 a21 = sin60*((a2+a5)-(a1+a4)) 1088 b11 = (b2+b5) + (b1+b4) 1089 b20 = (b(ia+i)+b3) - 0.5_wp*b11 1090 b21 = sin60*((b2+b5)-(b1+b4)) 1091 c(ja+j) = (a(ia+i)+a3) + a11 1092 d(ja+j) = (b(ia+i)+b3) + b11 1093 c(jc+j) = a20 - b21 1094 d(jc+j) = a21 + b20 1095 c(je+j) = a20 + b21 1096 d(je+j) = a21 - b20 1097 a11 = (a2-a5) + (a4-a1) 1098 a20 = (a(ia+i)-a3) - 0.5_wp*a11 1099 a21 = sin60*((a4-a1)-(a2-a5)) 1100 b11 = (b5-b2) - (b4-b1) 1101 b20 = (b3-b(ia+i)) - 0.5_wp*b11 1102 b21 = sin60*((b5-b2)+(b4-b1)) 1103 c(jb+j) = a20 - b21 1104 d(jb+j) = a21 - b20 1105 c(jd+j) = a11 + (a(ia+i)-a3) 1106 d(jd+j) = b11 + (b3-b(ia+i)) 1107 c(jf+j) = a20 + b21 1108 d(jf+j) = a21 + b20 1109 i = i + inc3 1110 j = j + inc4 1111 END DO 1112 ibase = ibase + inc1 1113 jbase = jbase + inc2 1114 END DO 1115 ibase = ibase + ijump 1116 ja = ja + jink 1117 jb = jb + jink 1118 jc = jc + jink 1119 jd = jd - jink 1120 je = je - jink 1121 jf = jf - jink 1122 END DO 1123 IF (jc>jd) GO TO 170 1124 140 CONTINUE 1125 jbase = 0 1126 DO l = 1, la 1127 i = ibase 1128 j = jbase 1129 !DIR$ IVDEP 1130 !CDIR NODEP 1131 !OCL NOVREC 1132 DO ijk = 1, lot 1133 c(ja+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) + sin60*(a(ib+i)-a(if+i)) 1134 d(ja+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) - sin60*(a(ic+i)+a(ie+i)) 1135 c(jb+j) = a(ia+i) - (a(ic+i)-a(ie+i)) 1136 d(jb+j) = a(id+i) - (a(ib+i)+a(if+i)) 1137 c(jc+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) - sin60*(a(ib+i)-a(if+i)) 1138 d(jc+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) + sin60*(a(ic+i)+a(ie+i)) 1139 i = i + inc3 1140 j = j + inc4 1141 END DO 1142 ibase = ibase + inc1 1143 jbase = jbase + inc2 1144 END DO 1145 1146 GO TO 170 1147 150 CONTINUE 1148 z = 1.0_wp/REAL(n,KIND=wp) 1149 zsin60 = z*sin60 1150 DO l = 1, la 1151 i = ibase 1152 j = jbase 1153 !DIR$ IVDEP 1154 !CDIR NODEP 1155 !OCL NOVREC 1156 DO ijk = 1, lot 1157 a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i)) 1158 c(ja+j) = z*((a(ia+i)+a(id+i))+a11) 1159 c(jc+j) = z*((a(ia+i)+a(id+i))-0.5_wp*a11) 1160 d(jc+j) = zsin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i))) 1161 a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i)) 1162 c(jb+j) = z*((a(ia+i)-a(id+i))-0.5_wp*a11) 1163 d(jb+j) = zsin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i))) 1164 c(jd+j) = z*((a(ia+i)-a(id+i))+a11) 1165 i = i + inc3 1166 j = j + inc4 1167 END DO 1168 ibase = ibase + inc1 1169 jbase = jbase + inc2 1170 END DO 1171 GO TO 170 1172 1173 ! Coding for factor 8 1174 1175 160 CONTINUE 1176 ibad = 3 1177 IF (la/=m) GO TO 180 1178 ia = 1 1179 ib = ia + iink 1180 ic = ib + iink 1181 id = ic + iink 1182 ie = id + iink 1183 if = ie + iink 1184 ig = if + iink 1185 ih = ig + iink 1186 ja = 1 1187 jb = ja + la*inc2 1188 jc = jb + 2*m*inc2 1189 jd = jc + 2*m*inc2 1190 je = jd + 2*m*inc2 1191 z = 1.0_wp/REAL(n,KIND=wp) 1192 zsin45 = z*SQRT(0.5_wp) 1193 1194 DO l = 1, la 1195 i = ibase 1196 j = jbase 1197 !DIR$ IVDEP 1198 !CDIR NODEP 1199 !OCL NOVREC 1200 DO ijk = 1, lot 1201 c(ja+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))+((a(id+i)+ & 1202 & a(ih+i))+(a(ib+i)+a(if+i)))) 1203 c(je+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))-((a(id+i)+ & 1204 & a(ih+i))+(a(ib+i)+a(if+i)))) 1205 c(jc+j) = z*((a(ia+i)+a(ie+i))-(a(ic+i)+a(ig+i))) 1206 d(jc+j) = z*((a(id+i)+a(ih+i))-(a(ib+i)+a(if+i))) 1207 c(jb+j) = z*(a(ia+i)-a(ie+i)) + zsin45*((a(ih+i)-a(id+i))-(a(if+ & 1208 & i)-a(ib+i))) 1209 c(jd+j) = z*(a(ia+i)-a(ie+i)) - zsin45*((a(ih+i)-a(id+i))-(a(if+ & 1210 & i)-a(ib+i))) 1211 d(jb+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) + & 1212 & z*(a(ig+i)-a(ic+i)) 1213 d(jd+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) - & 1214 & z*(a(ig+i)-a(ic+i)) 1215 i = i + inc3 1216 j = j + inc4 1217 END DO 1218 ibase = ibase + inc1 1219 jbase = jbase + inc2 1220 END DO 1221 1222 ! Return 1223 1224 170 CONTINUE 1225 ibad = 0 1226 180 CONTINUE 1227 ierr = ibad 1228 RETURN 1229 END SUBROUTINE qpassm 488 IF ( igo == 7 ) igo = 6 489 490 IF (igo < 1 .OR. igo > 6 ) THEN 491 ierr = 2 492 RETURN 493 ENDIF 494 495 496 SELECT CASE ( igo ) 497 498 ! 499 !-- Coding for factor 2 500 CASE ( 1 ) 501 ia = 1 502 ib = ia + iink 503 ja = 1 504 jb = ja + (2*m-la) * inc2 505 506 IF ( la /= m ) THEN 507 508 DO l = 1, la 509 i = ibase 510 j = jbase 511 !DIR$ IVDEP 512 !CDIR NODEP 513 !OCL NOVREC 514 DO ijk = 1, lot 515 c(ja+j) = a(ia+i) + a(ib+i) 516 c(jb+j) = a(ia+i) - a(ib+i) 517 i = i + inc3 518 j = j + inc4 519 ENDDO 520 ibase = ibase + inc1 521 jbase = jbase + inc2 522 ENDDO 523 524 ja = ja + jink 525 jink = 2 * jink 526 jb = jb - jink 527 ibase = ibase + ijump 528 ijump = 2 * ijump + iink 529 530 IF ( ja /= jb ) THEN 531 532 DO k = la, kstop, la 533 kb = k + k 534 c1 = trigs(kb+1) 535 s1 = trigs(kb+2) 536 jbase = 0 537 DO l = 1, la 538 i = ibase 539 j = jbase 540 !DIR$ IVDEP 541 !CDIR NODEP 542 !OCL NOVREC 543 DO ijk = 1, lot 544 c(ja+j) = a(ia+i) + (c1*a(ib+i)+s1*b(ib+i)) 545 c(jb+j) = a(ia+i) - (c1*a(ib+i)+s1*b(ib+i)) 546 d(ja+j) = (c1*b(ib+i)-s1*a(ib+i)) + b(ia+i) 547 d(jb+j) = (c1*b(ib+i)-s1*a(ib+i)) - b(ia+i) 548 i = i + inc3 549 j = j + inc4 550 ENDDO 551 ibase = ibase + inc1 552 jbase = jbase + inc2 553 ENDDO 554 555 ibase = ibase + ijump 556 ja = ja + jink 557 jb = jb - jink 558 ENDDO 559 560 IF ( ja > jb ) RETURN 561 562 ENDIF 563 564 jbase = 0 565 DO l = 1, la 566 i = ibase 567 j = jbase 568 !DIR$ IVDEP 569 !CDIR NODEP 570 !OCL NOVREC 571 DO ijk = 1, lot 572 c(ja+j) = a(ia+i) 573 d(ja+j) = -a(ib+i) 574 i = i + inc3 575 j = j + inc4 576 577 ENDDO 578 ibase = ibase + inc1 579 jbase = jbase + inc2 580 ENDDO 581 582 ELSE 583 584 z = 1.0_wp/REAL(n,KIND=wp) 585 DO l = 1, la 586 i = ibase 587 j = jbase 588 !DIR$ IVDEP 589 !CDIR NODEP 590 !OCL NOVREC 591 DO ijk = 1, lot 592 c(ja+j) = z*(a(ia+i)+a(ib+i)) 593 c(jb+j) = z*(a(ia+i)-a(ib+i)) 594 i = i + inc3 595 j = j + inc4 596 ENDDO 597 ibase = ibase + inc1 598 jbase = jbase + inc2 599 ENDDO 600 601 ENDIF 602 603 ! 604 !-- Coding for factor 3 605 CASE ( 2 ) 606 607 ia = 1 608 ib = ia + iink 609 ic = ib + iink 610 ja = 1 611 jb = ja + (2*m-la) * inc2 612 jc = jb 613 614 IF ( la /= m ) THEN 615 616 DO l = 1, la 617 i = ibase 618 j = jbase 619 !DIR$ IVDEP 620 !CDIR NODEP 621 !OCL NOVREC 622 DO ijk = 1, lot 623 c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i)) 624 c(jb+j) = a(ia+i) - 0.5_wp*(a(ib+i)+a(ic+i)) 625 d(jb+j) = sin60*(a(ic+i)-a(ib+i)) 626 i = i + inc3 627 j = j + inc4 628 ENDDO 629 ibase = ibase + inc1 630 jbase = jbase + inc2 631 ENDDO 632 633 ja = ja + jink 634 jink = 2 * jink 635 jb = jb + jink 636 jc = jc - jink 637 ibase = ibase + ijump 638 ijump = 2 * ijump + iink 639 640 IF ( ja /= jc ) THEN 641 642 DO k = la, kstop, la 643 kb = k + k 644 kc = kb + kb 645 c1 = trigs(kb+1) 646 s1 = trigs(kb+2) 647 c2 = trigs(kc+1) 648 s2 = trigs(kc+2) 649 650 jbase = 0 651 DO l = 1, la 652 i = ibase 653 j = jbase 654 !DIR$ IVDEP 655 !CDIR NODEP 656 !OCL NOVREC 657 DO ijk = 1, lot 658 a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c2*a(ic+i)+s2*b(ic+i)) 659 b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c2*b(ic+i)-s2*a(ic+i)) 660 a2 = a(ia+i) - 0.5_wp*a1 661 b2 = b(ia+i) - 0.5_wp*b1 662 a3 = sin60*((c1*a(ib+i)+s1*b(ib+i))-(c2*a(ic+i)+s2*b(ic+i))) 663 b3 = sin60*((c1*b(ib+i)-s1*a(ib+i))-(c2*b(ic+i)-s2*a(ic+i))) 664 c(ja+j) = a(ia+i) + a1 665 d(ja+j) = b(ia+i) + b1 666 c(jb+j) = a2 + b3 667 d(jb+j) = b2 - a3 668 c(jc+j) = a2 - b3 669 d(jc+j) = -(b2+a3) 670 i = i + inc3 671 j = j + inc4 672 ENDDO 673 ibase = ibase + inc1 674 jbase = jbase + inc2 675 ENDDO 676 ibase = ibase + ijump 677 ja = ja + jink 678 jb = jb + jink 679 jc = jc - jink 680 ENDDO 681 682 IF ( ja > jc ) RETURN 683 684 ENDIF 685 686 jbase = 0 687 DO l = 1, la 688 i = ibase 689 j = jbase 690 !DIR$ IVDEP 691 !CDIR NODEP 692 !OCL NOVREC 693 DO ijk = 1, lot 694 c(ja+j) = a(ia+i) + 0.5_wp*(a(ib+i)-a(ic+i)) 695 d(ja+j) = -sin60*(a(ib+i)+a(ic+i)) 696 c(jb+j) = a(ia+i) - (a(ib+i)-a(ic+i)) 697 i = i + inc3 698 j = j + inc4 699 ENDDO 700 ibase = ibase + inc1 701 jbase = jbase + inc2 702 ENDDO 703 704 ELSE 705 706 z = 1.0_wp / REAL( n, KIND=wp ) 707 zsin60 = z*sin60 708 DO l = 1, la 709 i = ibase 710 j = jbase 711 !DIR$ IVDEP 712 !CDIR NODEP 713 !OCL NOVREC 714 DO ijk = 1, lot 715 c(ja+j) = z*(a(ia+i)+(a(ib+i)+a(ic+i))) 716 c(jb+j) = z*(a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i))) 717 d(jb+j) = zsin60*(a(ic+i)-a(ib+i)) 718 i = i + inc3 719 j = j + inc4 720 ENDDO 721 ibase = ibase + inc1 722 jbase = jbase + inc2 723 ENDDO 724 725 ENDIF 726 727 ! 728 !-- Coding for factor 4 729 CASE ( 3 ) 730 731 ia = 1 732 ib = ia + iink 733 ic = ib + iink 734 id = ic + iink 735 ja = 1 736 jb = ja + (2*m-la) * inc2 737 jc = jb + 2*m*inc2 738 jd = jb 739 740 IF ( la /= m ) THEN 741 742 DO l = 1, la 743 744 i = ibase 745 j = jbase 746 !DIR$ IVDEP 747 !CDIR NODEP 748 !OCL NOVREC 749 DO ijk = 1, lot 750 c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i)) 751 c(jc+j) = (a(ia+i)+a(ic+i)) - (a(ib+i)+a(id+i)) 752 c(jb+j) = a(ia+i) - a(ic+i) 753 d(jb+j) = a(id+i) - a(ib+i) 754 i = i + inc3 755 j = j + inc4 756 ENDDO 757 ibase = ibase + inc1 758 jbase = jbase + inc2 759 ENDDO 760 761 ja = ja + jink 762 jink = 2 * jink 763 jb = jb + jink 764 jc = jc - jink 765 jd = jd - jink 766 ibase = ibase + ijump 767 ijump = 2 * ijump + iink 768 769 IF ( jb /= jc ) THEN 770 771 DO k = la, kstop, la 772 kb = k + k 773 kc = kb + kb 774 kd = kc + kb 775 c1 = trigs(kb+1) 776 s1 = trigs(kb+2) 777 c2 = trigs(kc+1) 778 s2 = trigs(kc+2) 779 c3 = trigs(kd+1) 780 s3 = trigs(kd+2) 781 jbase = 0 782 DO l = 1, la 783 i = ibase 784 j = jbase 785 !DIR$ IVDEP 786 !CDIR NODEP 787 !OCL NOVREC 788 DO ijk = 1, lot 789 a0 = a(ia+i) + (c2*a(ic+i)+s2*b(ic+i)) 790 a2 = a(ia+i) - (c2*a(ic+i)+s2*b(ic+i)) 791 a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c3*a(id+i)+s3*b(id+i)) 792 a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c3*a(id+i)+s3*b(id+i)) 793 b0 = b(ia+i) + (c2*b(ic+i)-s2*a(ic+i)) 794 b2 = b(ia+i) - (c2*b(ic+i)-s2*a(ic+i)) 795 b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c3*b(id+i)-s3*a(id+i)) 796 b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c3*b(id+i)-s3*a(id+i)) 797 c(ja+j) = a0 + a1 798 c(jc+j) = a0 - a1 799 d(ja+j) = b0 + b1 800 d(jc+j) = b1 - b0 801 c(jb+j) = a2 + b3 802 c(jd+j) = a2 - b3 803 d(jb+j) = b2 - a3 804 d(jd+j) = -(b2+a3) 805 i = i + inc3 806 j = j + inc4 807 ENDDO 808 ibase = ibase + inc1 809 jbase = jbase + inc2 810 ENDDO 811 ibase = ibase + ijump 812 ja = ja + jink 813 jb = jb + jink 814 jc = jc - jink 815 jd = jd - jink 816 ENDDO 817 818 IF ( jb > jc ) RETURN 819 820 ENDIF 821 822 sin45 = SQRT( 0.5_wp ) 823 jbase = 0 824 DO l = 1, la 825 i = ibase 826 j = jbase 827 !DIR$ IVDEP 828 !CDIR NODEP 829 !OCL NOVREC 830 DO ijk = 1, lot 831 c(ja+j) = a(ia+i) + sin45*(a(ib+i)-a(id+i)) 832 c(jb+j) = a(ia+i) - sin45*(a(ib+i)-a(id+i)) 833 d(ja+j) = -a(ic+i) - sin45*(a(ib+i)+a(id+i)) 834 d(jb+j) = a(ic+i) - sin45*(a(ib+i)+a(id+i)) 835 i = i + inc3 836 j = j + inc4 837 ENDDO 838 ibase = ibase + inc1 839 jbase = jbase + inc2 840 ENDDO 841 842 ELSE 843 844 z = 1.0_wp / REAL( n, KIND=wp ) 845 DO l = 1, la 846 i = ibase 847 j = jbase 848 !DIR$ IVDEP 849 !CDIR NODEP 850 !OCL NOVREC 851 DO ijk = 1, lot 852 c(ja+j) = z*((a(ia+i)+a(ic+i))+(a(ib+i)+a(id+i))) 853 c(jc+j) = z*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) 854 c(jb+j) = z*(a(ia+i)-a(ic+i)) 855 d(jb+j) = z*(a(id+i)-a(ib+i)) 856 i = i + inc3 857 j = j + inc4 858 ENDDO 859 ibase = ibase + inc1 860 jbase = jbase + inc2 861 ENDDO 862 863 ENDIF 864 865 ! 866 !-- Coding for factor 5 867 CASE ( 4 ) 868 869 ia = 1 870 ib = ia + iink 871 ic = ib + iink 872 id = ic + iink 873 ie = id + iink 874 ja = 1 875 jb = ja + (2*m-la) * inc2 876 jc = jb + 2*m*inc2 877 jd = jc 878 je = jb 879 880 IF ( la /= m ) THEN 881 882 DO l = 1, la 883 i = ibase 884 j = jbase 885 !DIR$ IVDEP 886 !CDIR NODEP 887 !OCL NOVREC 888 DO ijk = 1, lot 889 a1 = a(ib+i) + a(ie+i) 890 a3 = a(ib+i) - a(ie+i) 891 a2 = a(ic+i) + a(id+i) 892 a4 = a(ic+i) - a(id+i) 893 a5 = a(ia+i) - 0.25_wp*(a1+a2) 894 a6 = qrt5*(a1-a2) 895 c(ja+j) = a(ia+i) + (a1+a2) 896 c(jb+j) = a5 + a6 897 c(jc+j) = a5 - a6 898 d(jb+j) = -sin72*a3 - sin36*a4 899 d(jc+j) = -sin36*a3 + sin72*a4 900 i = i + inc3 901 j = j + inc4 902 ENDDO 903 ibase = ibase + inc1 904 jbase = jbase + inc2 905 ENDDO 906 907 ja = ja + jink 908 jink = 2 * jink 909 jb = jb + jink 910 jc = jc + jink 911 jd = jd - jink 912 je = je - jink 913 ibase = ibase + ijump 914 ijump = 2 * ijump + iink 915 916 IF ( jb /= jd ) THEN 917 918 DO k = la, kstop, la 919 kb = k + k 920 kc = kb + kb 921 kd = kc + kb 922 ke = kd + kb 923 c1 = trigs(kb+1) 924 s1 = trigs(kb+2) 925 c2 = trigs(kc+1) 926 s2 = trigs(kc+2) 927 c3 = trigs(kd+1) 928 s3 = trigs(kd+2) 929 c4 = trigs(ke+1) 930 s4 = trigs(ke+2) 931 jbase = 0 932 DO l = 1, la 933 i = ibase 934 j = jbase 935 !DIR$ IVDEP 936 !CDIR NODEP 937 !OCL NOVREC 938 DO ijk = 1, lot 939 a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c4*a(ie+i)+s4*b(ie+i)) 940 a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c4*a(ie+i)+s4*b(ie+i)) 941 a2 = (c2*a(ic+i)+s2*b(ic+i)) + (c3*a(id+i)+s3*b(id+i)) 942 a4 = (c2*a(ic+i)+s2*b(ic+i)) - (c3*a(id+i)+s3*b(id+i)) 943 b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c4*b(ie+i)-s4*a(ie+i)) 944 b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c4*b(ie+i)-s4*a(ie+i)) 945 b2 = (c2*b(ic+i)-s2*a(ic+i)) + (c3*b(id+i)-s3*a(id+i)) 946 b4 = (c2*b(ic+i)-s2*a(ic+i)) - (c3*b(id+i)-s3*a(id+i)) 947 a5 = a(ia+i) - 0.25_wp*(a1+a2) 948 a6 = qrt5*(a1-a2) 949 b5 = b(ia+i) - 0.25_wp*(b1+b2) 950 b6 = qrt5*(b1-b2) 951 a10 = a5 + a6 952 a20 = a5 - a6 953 b10 = b5 + b6 954 b20 = b5 - b6 955 a11 = sin72*b3 + sin36*b4 956 a21 = sin36*b3 - sin72*b4 957 b11 = sin72*a3 + sin36*a4 958 b21 = sin36*a3 - sin72*a4 959 c(ja+j) = a(ia+i) + (a1+a2) 960 c(jb+j) = a10 + a11 961 c(je+j) = a10 - a11 962 c(jc+j) = a20 + a21 963 c(jd+j) = a20 - a21 964 d(ja+j) = b(ia+i) + (b1+b2) 965 d(jb+j) = b10 - b11 966 d(je+j) = -(b10+b11) 967 d(jc+j) = b20 - b21 968 d(jd+j) = -(b20+b21) 969 i = i + inc3 970 j = j + inc4 971 ENDDO 972 ibase = ibase + inc1 973 jbase = jbase + inc2 974 ENDDO 975 ibase = ibase + ijump 976 ja = ja + jink 977 jb = jb + jink 978 jc = jc + jink 979 jd = jd - jink 980 je = je - jink 981 ENDDO 982 983 IF ( jb > jd ) RETURN 984 985 ENDIF 986 987 jbase = 0 988 DO l = 1, la 989 i = ibase 990 j = jbase 991 !DIR$ IVDEP 992 !CDIR NODEP 993 !OCL NOVREC 994 DO ijk = 1, lot 995 a1 = a(ib+i) + a(ie+i) 996 a3 = a(ib+i) - a(ie+i) 997 a2 = a(ic+i) + a(id+i) 998 a4 = a(ic+i) - a(id+i) 999 a5 = a(ia+i) + 0.25_wp*(a3-a4) 1000 a6 = qrt5*(a3+a4) 1001 c(ja+j) = a5 + a6 1002 c(jb+j) = a5 - a6 1003 c(jc+j) = a(ia+i) - (a3-a4) 1004 d(ja+j) = -sin36*a1 - sin72*a2 1005 d(jb+j) = -sin72*a1 + sin36*a2 1006 i = i + inc3 1007 j = j + inc4 1008 ENDDO 1009 ibase = ibase + inc1 1010 jbase = jbase + inc2 1011 ENDDO 1012 1013 ELSE 1014 1015 z = 1.0_wp / REAL( n, KIND=wp ) 1016 zqrt5 = z * qrt5 1017 zsin36 = z * sin36 1018 zsin72 = z * sin72 1019 DO l = 1, la 1020 i = ibase 1021 j = jbase 1022 !DIR$ IVDEP 1023 !CDIR NODEP 1024 !OCL NOVREC 1025 DO ijk = 1, lot 1026 a1 = a(ib+i) + a(ie+i) 1027 a3 = a(ib+i) - a(ie+i) 1028 a2 = a(ic+i) + a(id+i) 1029 a4 = a(ic+i) - a(id+i) 1030 a5 = z*(a(ia+i)-0.25_wp*(a1+a2)) 1031 a6 = zqrt5*(a1-a2) 1032 c(ja+j) = z*(a(ia+i)+(a1+a2)) 1033 c(jb+j) = a5 + a6 1034 c(jc+j) = a5 - a6 1035 d(jb+j) = -zsin72*a3 - zsin36*a4 1036 d(jc+j) = -zsin36*a3 + zsin72*a4 1037 i = i + inc3 1038 j = j + inc4 1039 ENDDO 1040 ibase = ibase + inc1 1041 jbase = jbase + inc2 1042 ENDDO 1043 1044 ENDIF 1045 1046 ! 1047 !-- Coding for factor 6 1048 CASE ( 5 ) 1049 1050 ia = 1 1051 ib = ia + iink 1052 ic = ib + iink 1053 id = ic + iink 1054 ie = id + iink 1055 if = ie + iink 1056 ja = 1 1057 jb = ja + (2*m-la) * inc2 1058 jc = jb + 2*m*inc2 1059 jd = jc + 2*m*inc2 1060 je = jc 1061 jf = jb 1062 1063 IF ( la /= m ) THEN 1064 1065 DO l = 1, la 1066 i = ibase 1067 j = jbase 1068 !DIR$ IVDEP 1069 !CDIR NODEP 1070 !OCL NOVREC 1071 DO ijk = 1, lot 1072 a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i)) 1073 c(ja+j) = (a(ia+i)+a(id+i)) + a11 1074 c(jc+j) = (a(ia+i)+a(id+i)-0.5_wp*a11) 1075 d(jc+j) = sin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i))) 1076 a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i)) 1077 c(jb+j) = (a(ia+i)-a(id+i)) - 0.5_wp*a11 1078 d(jb+j) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i))) 1079 c(jd+j) = (a(ia+i)-a(id+i)) + a11 1080 i = i + inc3 1081 j = j + inc4 1082 ENDDO 1083 ibase = ibase + inc1 1084 jbase = jbase + inc2 1085 ENDDO 1086 ja = ja + jink 1087 jink = 2 * jink 1088 jb = jb + jink 1089 jc = jc + jink 1090 jd = jd - jink 1091 je = je - jink 1092 jf = jf - jink 1093 ibase = ibase + ijump 1094 ijump = 2 * ijump + iink 1095 1096 IF ( jc /= jd ) THEN 1097 1098 DO k = la, kstop, la 1099 kb = k + k 1100 kc = kb + kb 1101 kd = kc + kb 1102 ke = kd + kb 1103 kf = ke + kb 1104 c1 = trigs(kb+1) 1105 s1 = trigs(kb+2) 1106 c2 = trigs(kc+1) 1107 s2 = trigs(kc+2) 1108 c3 = trigs(kd+1) 1109 s3 = trigs(kd+2) 1110 c4 = trigs(ke+1) 1111 s4 = trigs(ke+2) 1112 c5 = trigs(kf+1) 1113 s5 = trigs(kf+2) 1114 jbase = 0 1115 DO l = 1, la 1116 i = ibase 1117 j = jbase 1118 !DIR$ IVDEP 1119 !CDIR NODEP 1120 !OCL NOVREC 1121 DO ijk = 1, lot 1122 a1 = c1*a(ib+i) + s1*b(ib+i) 1123 b1 = c1*b(ib+i) - s1*a(ib+i) 1124 a2 = c2*a(ic+i) + s2*b(ic+i) 1125 b2 = c2*b(ic+i) - s2*a(ic+i) 1126 a3 = c3*a(id+i) + s3*b(id+i) 1127 b3 = c3*b(id+i) - s3*a(id+i) 1128 a4 = c4*a(ie+i) + s4*b(ie+i) 1129 b4 = c4*b(ie+i) - s4*a(ie+i) 1130 a5 = c5*a(if+i) + s5*b(if+i) 1131 b5 = c5*b(if+i) - s5*a(if+i) 1132 a11 = (a2+a5) + (a1+a4) 1133 a20 = (a(ia+i)+a3) - 0.5_wp*a11 1134 a21 = sin60*((a2+a5)-(a1+a4)) 1135 b11 = (b2+b5) + (b1+b4) 1136 b20 = (b(ia+i)+b3) - 0.5_wp*b11 1137 b21 = sin60*((b2+b5)-(b1+b4)) 1138 c(ja+j) = (a(ia+i)+a3) + a11 1139 d(ja+j) = (b(ia+i)+b3) + b11 1140 c(jc+j) = a20 - b21 1141 d(jc+j) = a21 + b20 1142 c(je+j) = a20 + b21 1143 d(je+j) = a21 - b20 1144 a11 = (a2-a5) + (a4-a1) 1145 a20 = (a(ia+i)-a3) - 0.5_wp*a11 1146 a21 = sin60*((a4-a1)-(a2-a5)) 1147 b11 = (b5-b2) - (b4-b1) 1148 b20 = (b3-b(ia+i)) - 0.5_wp*b11 1149 b21 = sin60*((b5-b2)+(b4-b1)) 1150 c(jb+j) = a20 - b21 1151 d(jb+j) = a21 - b20 1152 c(jd+j) = a11 + (a(ia+i)-a3) 1153 d(jd+j) = b11 + (b3-b(ia+i)) 1154 c(jf+j) = a20 + b21 1155 d(jf+j) = a21 + b20 1156 i = i + inc3 1157 j = j + inc4 1158 ENDDO 1159 ibase = ibase + inc1 1160 jbase = jbase + inc2 1161 ENDDO 1162 ibase = ibase + ijump 1163 ja = ja + jink 1164 jb = jb + jink 1165 jc = jc + jink 1166 jd = jd - jink 1167 je = je - jink 1168 jf = jf - jink 1169 ENDDO 1170 1171 IF ( jc > jd ) RETURN 1172 1173 ENDIF 1174 1175 jbase = 0 1176 DO l = 1, la 1177 i = ibase 1178 j = jbase 1179 !DIR$ IVDEP 1180 !CDIR NODEP 1181 !OCL NOVREC 1182 DO ijk = 1, lot 1183 c(ja+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) + sin60*(a(ib+i)-a(if+i)) 1184 d(ja+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) - sin60*(a(ic+i)+a(ie+i)) 1185 c(jb+j) = a(ia+i) - (a(ic+i)-a(ie+i)) 1186 d(jb+j) = a(id+i) - (a(ib+i)+a(if+i)) 1187 c(jc+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) - sin60*(a(ib+i)-a(if+i)) 1188 d(jc+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) + sin60*(a(ic+i)+a(ie+i)) 1189 i = i + inc3 1190 j = j + inc4 1191 ENDDO 1192 ibase = ibase + inc1 1193 jbase = jbase + inc2 1194 ENDDO 1195 1196 ELSE 1197 1198 z = 1.0_wp/REAL(n,KIND=wp) 1199 zsin60 = z*sin60 1200 DO l = 1, la 1201 i = ibase 1202 j = jbase 1203 !DIR$ IVDEP 1204 !CDIR NODEP 1205 !OCL NOVREC 1206 DO ijk = 1, lot 1207 a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i)) 1208 c(ja+j) = z*((a(ia+i)+a(id+i))+a11) 1209 c(jc+j) = z*((a(ia+i)+a(id+i))-0.5_wp*a11) 1210 d(jc+j) = zsin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i))) 1211 a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i)) 1212 c(jb+j) = z*((a(ia+i)-a(id+i))-0.5_wp*a11) 1213 d(jb+j) = zsin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i))) 1214 c(jd+j) = z*((a(ia+i)-a(id+i))+a11) 1215 i = i + inc3 1216 j = j + inc4 1217 ENDDO 1218 ibase = ibase + inc1 1219 jbase = jbase + inc2 1220 ENDDO 1221 1222 ENDIF 1223 1224 ! 1225 !-- Coding for factor 8 1226 CASE ( 6 ) 1227 1228 IF ( la /= m ) THEN 1229 ierr = 3 1230 RETURN 1231 ENDIF 1232 1233 ia = 1 1234 ib = ia + iink 1235 ic = ib + iink 1236 id = ic + iink 1237 ie = id + iink 1238 if = ie + iink 1239 ig = if + iink 1240 ih = ig + iink 1241 ja = 1 1242 jb = ja + la * inc2 1243 jc = jb + 2*m*inc2 1244 jd = jc + 2*m*inc2 1245 je = jd + 2*m*inc2 1246 z = 1.0_wp / REAL( n, KIND=wp ) 1247 zsin45 = z * SQRT( 0.5_wp ) 1248 1249 DO l = 1, la 1250 i = ibase 1251 j = jbase 1252 1253 !DIR$ IVDEP 1254 !CDIR NODEP 1255 !OCL NOVREC 1256 DO ijk = 1, lot 1257 c(ja+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))+((a(id+i)+ a(ih+i))+(a(ib+i)+a(if+i)))) 1258 c(je+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))-((a(id+i)+ a(ih+i))+(a(ib+i)+a(if+i)))) 1259 c(jc+j) = z*((a(ia+i)+a(ie+i))-(a(ic+i)+a(ig+i))) 1260 d(jc+j) = z*((a(id+i)+a(ih+i))-(a(ib+i)+a(if+i))) 1261 c(jb+j) = z*(a(ia+i)-a(ie+i)) + zsin45*((a(ih+i)-a(id+i))-(a(if+i)-a(ib+i))) 1262 c(jd+j) = z*(a(ia+i)-a(ie+i)) - zsin45*((a(ih+i)-a(id+i))-(a(if+i)-a(ib+i))) 1263 d(jb+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) + z*(a(ig+i)-a(ic+i)) 1264 d(jd+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) - z*(a(ig+i)-a(ic+i)) 1265 i = i + inc3 1266 j = j + inc4 1267 ENDDO 1268 ibase = ibase + inc1 1269 jbase = jbase + inc2 1270 ENDDO 1271 1272 END SELECT 1273 1274 END SUBROUTINE qpassm 1230 1275 1231 1276 !------------------------------------------------------------------------------! 1232 1277 ! Description: 1233 1278 ! ------------ 1234 !> @todo Missing subroutine description.1279 !> Same as qpassm, but for backward fft 1235 1280 !------------------------------------------------------------------------------! 1236 SUBROUTINE rpassm(a,b,c,d,trigs,inc1,inc2,inc3,inc4,lot,n,ifac,la,ierr) 1237 ! Dimension a(n),b(n),c(n),d(n),trigs(n) 1281 SUBROUTINE rpassm( a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr ) 1238 1282 1239 1283 USE kinds … … 1241 1285 IMPLICIT NONE 1242 1286 1243 ! Scalar arguments1244 1287 INTEGER(iwp) :: ierr !< 1245 1288 INTEGER(iwp) :: ifac !< … … 1251 1294 INTEGER(iwp) :: lot !< 1252 1295 INTEGER(iwp) :: n !< 1253 1254 ! Array arguments 1296 ! 1297 !-- Arrays are dimensioned with n 1255 1298 REAL(wp) :: a(*) !< 1256 1299 REAL(wp) :: b(*) !< … … 1259 1302 REAL(wp) :: trigs(*) !< 1260 1303 1261 ! Local scalars:1262 1304 REAL(wp) :: c1 !< 1263 1305 REAL(wp) :: c2 !< … … 1284 1326 INTEGER(iwp) :: ia !< 1285 1327 INTEGER(iwp) :: ib !< 1286 INTEGER(iwp) :: ibad !<1287 1328 INTEGER(iwp) :: ibase !< 1288 1329 INTEGER(iwp) :: ic !< … … 1315 1356 INTEGER(iwp) :: m !< 1316 1357 1317 ! Local arrays:1318 1358 REAL(wp) :: a10(nfft) !< 1319 1359 REAL(wp) :: a11(nfft) !< … … 1325 1365 REAL(wp) :: b21(nfft) !< 1326 1366 1327 ! Intrinsic functions 1328 ! INTRINSIC SQRT 1329 1330 ! Data statements 1331 DATA sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/, & 1332 & qrt5/0.559016994374947_wp/, sin60/0.866025403784437_wp/ 1333 1334 1335 ! Executable statements 1336 1337 m = n/ifac 1338 iink = la*inc1 1339 jink = la*inc2 1340 jump = (ifac-1)*jink 1341 kstop = (n-ifac)/(2*ifac) 1342 1343 ibad = 1 1344 IF (lot>nfft) GO TO 180 1367 1368 DATA sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/, & 1369 qrt5/0.559016994374947_wp/, sin60/0.866025403784437_wp/ 1370 1371 1372 ierr = 0 1373 1374 m = n / ifac 1375 iink = la * inc1 1376 jink = la * inc2 1377 jump = (ifac-1) * jink 1378 kstop = (n-ifac) / (2*ifac) 1379 1380 IF ( lot > nfft ) THEN 1381 ierr = 1 1382 RETURN 1383 ENDIF 1345 1384 ibase = 0 1346 1385 jbase = 0 1347 1386 igo = ifac - 1 1348 IF (igo==7) igo = 6 1349 ibad = 2 1350 IF (igo<1 .OR. igo>6) GO TO 180 1351 GO TO (10,40,70,100,130,160) igo 1352 1353 ! Coding for factor 2 1354 1355 10 CONTINUE 1356 ia = 1 1357 ib = ia + (2*m-la)*inc1 1358 ja = 1 1359 jb = ja + jink 1360 1361 IF (la==m) GO TO 30 1362 1363 DO l = 1, la 1364 i = ibase 1365 j = jbase 1366 !DIR$ IVDEP 1367 !CDIR NODEP 1368 !OCL NOVREC 1369 DO ijk = 1, lot 1370 c(ja+j) = a(ia+i) + a(ib+i) 1371 c(jb+j) = a(ia+i) - a(ib+i) 1372 i = i + inc3 1373 j = j + inc4 1374 END DO 1375 ibase = ibase + inc1 1376 jbase = jbase + inc2 1377 END DO 1378 ia = ia + iink 1379 iink = 2*iink 1380 ib = ib - iink 1381 ibase = 0 1382 jbase = jbase + jump 1383 jump = 2*jump + jink 1384 IF (ia==ib) GO TO 20 1385 DO k = la, kstop, la 1386 kb = k + k 1387 c1 = trigs(kb+1) 1388 s1 = trigs(kb+2) 1389 ibase = 0 1390 DO l = 1, la 1391 i = ibase 1392 j = jbase 1393 !DIR$ IVDEP 1394 !CDIR NODEP 1395 !OCL NOVREC 1396 DO ijk = 1, lot 1397 c(ja+j) = a(ia+i) + a(ib+i) 1398 d(ja+j) = b(ia+i) - b(ib+i) 1399 c(jb+j) = c1*(a(ia+i)-a(ib+i)) - s1*(b(ia+i)+b(ib+i)) 1400 d(jb+j) = s1*(a(ia+i)-a(ib+i)) + c1*(b(ia+i)+b(ib+i)) 1401 i = i + inc3 1402 j = j + inc4 1403 END DO 1404 ibase = ibase + inc1 1405 jbase = jbase + inc2 1406 END DO 1407 ia = ia + iink 1408 ib = ib - iink 1409 jbase = jbase + jump 1410 END DO 1411 IF (ia>ib) GO TO 170 1412 20 CONTINUE 1413 ibase = 0 1414 DO l = 1, la 1415 i = ibase 1416 j = jbase 1417 !DIR$ IVDEP 1418 !CDIR NODEP 1419 !OCL NOVREC 1420 DO ijk = 1, lot 1421 c(ja+j) = a(ia+i) 1422 c(jb+j) = -b(ia+i) 1423 i = i + inc3 1424 j = j + inc4 1425 END DO 1426 ibase = ibase + inc1 1427 jbase = jbase + inc2 1428 END DO 1429 1430 GO TO 170 1431 30 CONTINUE 1432 DO l = 1, la 1433 i = ibase 1434 j = jbase 1435 !DIR$ IVDEP 1436 !CDIR NODEP 1437 !OCL NOVREC 1438 DO ijk = 1, lot 1439 c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i)) 1440 c(jb+j) = 2.0_wp*(a(ia+i)-a(ib+i)) 1441 i = i + inc3 1442 j = j + inc4 1443 END DO 1444 ibase = ibase + inc1 1445 jbase = jbase + inc2 1446 END DO 1447 GO TO 170 1448 1449 ! Coding for factor 3 1450 1451 40 CONTINUE 1452 ia = 1 1453 ib = ia + (2*m-la)*inc1 1454 ic = ib 1455 ja = 1 1456 jb = ja + jink 1457 jc = jb + jink 1458 1459 IF (la==m) GO TO 60 1460 1461 DO l = 1, la 1462 i = ibase 1463 j = jbase 1464 !DIR$ IVDEP 1465 !CDIR NODEP 1466 !OCL NOVREC 1467 DO ijk = 1, lot 1468 c(ja+j) = a(ia+i) + a(ib+i) 1469 c(jb+j) = (a(ia+i)-0.5_wp*a(ib+i)) - (sin60*(b(ib+i))) 1470 c(jc+j) = (a(ia+i)-0.5_wp*a(ib+i)) + (sin60*(b(ib+i))) 1471 i = i + inc3 1472 j = j + inc4 1473 END DO 1474 ibase = ibase + inc1 1475 jbase = jbase + inc2 1476 END DO 1477 ia = ia + iink 1478 iink = 2*iink 1479 ib = ib + iink 1480 ic = ic - iink 1481 jbase = jbase + jump 1482 jump = 2*jump + jink 1483 IF (ia==ic) GO TO 50 1484 DO k = la, kstop, la 1485 kb = k + k 1486 kc = kb + kb 1487 c1 = trigs(kb+1) 1488 s1 = trigs(kb+2) 1489 c2 = trigs(kc+1) 1490 s2 = trigs(kc+2) 1491 ibase = 0 1492 DO l = 1, la 1493 i = ibase 1494 j = jbase 1495 !DIR$ IVDEP 1496 !CDIR NODEP 1497 !OCL NOVREC 1498 DO ijk = 1, lot 1499 c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i)) 1500 d(ja+j) = b(ia+i) + (b(ib+i)-b(ic+i)) 1501 c(jb+j) = c1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ & 1502 & b(ic+i)))) & 1503 & - s1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- & 1504 & a(ic+i)))) 1505 d(jb+j) = s1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ & 1506 & b(ic+i)))) & 1507 & + c1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- & 1508 & a(ic+i)))) 1509 c(jc+j) = c2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ & 1510 & b(ic+i)))) & 1511 & - s2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- & 1512 & a(ic+i)))) 1513 d(jc+j) = s2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ & 1514 & b(ic+i)))) & 1515 & + c2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- & 1516 & a(ic+i)))) 1517 i = i + inc3 1518 j = j + inc4 1519 END DO 1520 ibase = ibase + inc1 1521 jbase = jbase + inc2 1522 END DO 1523 ia = ia + iink 1524 ib = ib + iink 1525 ic = ic - iink 1526 jbase = jbase + jump 1527 END DO 1528 IF (ia>ic) GO TO 170 1529 50 CONTINUE 1530 ibase = 0 1531 DO l = 1, la 1532 i = ibase 1533 j = jbase 1534 !DIR$ IVDEP 1535 !CDIR NODEP 1536 !OCL NOVREC 1537 DO ijk = 1, lot 1538 c(ja+j) = a(ia+i) + a(ib+i) 1539 c(jb+j) = (0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i)) 1540 c(jc+j) = -(0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i)) 1541 i = i + inc3 1542 j = j + inc4 1543 END DO 1544 ibase = ibase + inc1 1545 jbase = jbase + inc2 1546 END DO 1547 1548 GO TO 170 1549 60 CONTINUE 1550 ssin60 = 2.0_wp*sin60 1551 DO l = 1, la 1552 i = ibase 1553 j = jbase 1554 !DIR$ IVDEP 1555 !CDIR NODEP 1556 !OCL NOVREC 1557 DO ijk = 1, lot 1558 c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i)) 1559 c(jb+j) = (2.0_wp*a(ia+i)-a(ib+i)) - (ssin60*b(ib+i)) 1560 c(jc+j) = (2.0_wp*a(ia+i)-a(ib+i)) + (ssin60*b(ib+i)) 1561 i = i + inc3 1562 j = j + inc4 1563 END DO 1564 ibase = ibase + inc1 1565 jbase = jbase + inc2 1566 END DO 1567 GO TO 170 1568 1569 ! Coding for factor 4 1570 1571 70 CONTINUE 1572 ia = 1 1573 ib = ia + (2*m-la)*inc1 1574 ic = ib + 2*m*inc1 1575 id = ib 1576 ja = 1 1577 jb = ja + jink 1578 jc = jb + jink 1579 jd = jc + jink 1580 1581 IF (la==m) GO TO 90 1582 1583 DO l = 1, la 1584 i = ibase 1585 j = jbase 1586 !DIR$ IVDEP 1587 !CDIR NODEP 1588 !OCL NOVREC 1589 DO ijk = 1, lot 1590 c(ja+j) = (a(ia+i)+a(ic+i)) + a(ib+i) 1591 c(jb+j) = (a(ia+i)-a(ic+i)) - b(ib+i) 1592 c(jc+j) = (a(ia+i)+a(ic+i)) - a(ib+i) 1593 c(jd+j) = (a(ia+i)-a(ic+i)) + b(ib+i) 1594 i = i + inc3 1595 j = j + inc4 1596 END DO 1597 ibase = ibase + inc1 1598 jbase = jbase + inc2 1599 END DO 1600 ia = ia + iink 1601 iink = 2*iink 1602 ib = ib + iink 1603 ic = ic - iink 1604 id = id - iink 1605 jbase = jbase + jump 1606 jump = 2*jump + jink 1607 IF (ib==ic) GO TO 80 1608 DO k = la, kstop, la 1609 kb = k + k 1610 kc = kb + kb 1611 kd = kc + kb 1612 c1 = trigs(kb+1) 1613 s1 = trigs(kb+2) 1614 c2 = trigs(kc+1) 1615 s2 = trigs(kc+2) 1616 c3 = trigs(kd+1) 1617 s3 = trigs(kd+2) 1618 ibase = 0 1619 DO l = 1, la 1620 i = ibase 1621 j = jbase 1622 !DIR$ IVDEP 1623 !CDIR NODEP 1624 !OCL NOVREC 1625 DO ijk = 1, lot 1626 c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i)) 1627 d(ja+j) = (b(ia+i)-b(ic+i)) + (b(ib+i)-b(id+i)) 1628 c(jc+j) = c2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) - s2*((b(ia+ & 1629 & i)-b(ic+i))-(b(ib+i)-b(id+i))) 1630 d(jc+j) = s2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) + c2*((b(ia+ & 1631 & i)-b(ic+i))-(b(ib+i)-b(id+i))) 1632 c(jb+j) = c1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) - s1*((b(ia+ & 1633 & i)+b(ic+i))+(a(ib+i)-a(id+i))) 1634 d(jb+j) = s1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) + c1*((b(ia+ & 1635 & i)+b(ic+i))+(a(ib+i)-a(id+i))) 1636 c(jd+j) = c3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) - s3*((b(ia+ & 1637 & i)+b(ic+i))-(a(ib+i)-a(id+i))) 1638 d(jd+j) = s3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) + c3*((b(ia+ & 1639 & i)+b(ic+i))-(a(ib+i)-a(id+i))) 1640 i = i + inc3 1641 j = j + inc4 1642 END DO 1643 ibase = ibase + inc1 1644 jbase = jbase + inc2 1645 END DO 1646 ia = ia + iink 1647 ib = ib + iink 1648 ic = ic - iink 1649 id = id - iink 1650 jbase = jbase + jump 1651 END DO 1652 IF (ib>ic) GO TO 170 1653 80 CONTINUE 1654 ibase = 0 1655 sin45 = SQRT(0.5_wp) 1656 DO l = 1, la 1657 i = ibase 1658 j = jbase 1659 !DIR$ IVDEP 1660 !CDIR NODEP 1661 !OCL NOVREC 1662 DO ijk = 1, lot 1663 c(ja+j) = a(ia+i) + a(ib+i) 1664 c(jb+j) = sin45*((a(ia+i)-a(ib+i))-(b(ia+i)+b(ib+i))) 1665 c(jc+j) = b(ib+i) - b(ia+i) 1666 c(jd+j) = -sin45*((a(ia+i)-a(ib+i))+(b(ia+i)+b(ib+i))) 1667 i = i + inc3 1668 j = j + inc4 1669 END DO 1670 ibase = ibase + inc1 1671 jbase = jbase + inc2 1672 END DO 1673 1674 GO TO 170 1675 90 CONTINUE 1676 DO l = 1, la 1677 i = ibase 1678 j = jbase 1679 !DIR$ IVDEP 1680 !CDIR NODEP 1681 !OCL NOVREC 1682 DO ijk = 1, lot 1683 c(ja+j) = 2.0_wp*((a(ia+i)+a(ic+i))+a(ib+i)) 1684 c(jb+j) = 2.0_wp*((a(ia+i)-a(ic+i))-b(ib+i)) 1685 c(jc+j) = 2.0_wp*((a(ia+i)+a(ic+i))-a(ib+i)) 1686 c(jd+j) = 2.0_wp*((a(ia+i)-a(ic+i))+b(ib+i)) 1687 i = i + inc3 1688 j = j + inc4 1689 END DO 1690 ibase = ibase + inc1 1691 jbase = jbase + inc2 1692 END DO 1693 1694 ! Coding for factor 5 1695 1696 GO TO 170 1697 100 CONTINUE 1698 ia = 1 1699 ib = ia + (2*m-la)*inc1 1700 ic = ib + 2*m*inc1 1701 id = ic 1702 ie = ib 1703 ja = 1 1704 jb = ja + jink 1705 jc = jb + jink 1706 jd = jc + jink 1707 je = jd + jink 1708 1709 IF (la==m) GO TO 120 1710 1711 DO l = 1, la 1712 i = ibase 1713 j = jbase 1714 !DIR$ IVDEP 1715 !CDIR NODEP 1716 !OCL NOVREC 1717 DO ijk = 1, lot 1718 c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i)) 1719 c(jb+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) - & 1720 & (sin72*b(ib+i)+sin36*b(ic+i)) 1721 c(jc+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) - & 1722 & (sin36*b(ib+i)-sin72*b(ic+i)) 1723 c(jd+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) + & 1724 & (sin36*b(ib+i)-sin72*b(ic+i)) 1725 c(je+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) + & 1726 & (sin72*b(ib+i)+sin36*b(ic+i)) 1727 i = i + inc3 1728 j = j + inc4 1729 END DO 1730 ibase = ibase + inc1 1731 jbase = jbase + inc2 1732 END DO 1733 ia = ia + iink 1734 iink = 2*iink 1735 ib = ib + iink 1736 ic = ic + iink 1737 id = id - iink 1738 ie = ie - iink 1739 jbase = jbase + jump 1740 jump = 2*jump + jink 1741 IF (ib==id) GO TO 110 1742 DO k = la, kstop, la 1743 kb = k + k 1744 kc = kb + kb 1745 kd = kc + kb 1746 ke = kd + kb 1747 c1 = trigs(kb+1) 1748 s1 = trigs(kb+2) 1749 c2 = trigs(kc+1) 1750 s2 = trigs(kc+2) 1751 c3 = trigs(kd+1) 1752 s3 = trigs(kd+2) 1753 c4 = trigs(ke+1) 1754 s4 = trigs(ke+2) 1755 ibase = 0 1756 DO l = 1, la 1757 i = ibase 1758 j = jbase 1759 !DIR$ IVDEP 1760 !CDIR NODEP 1761 !OCL NOVREC 1762 DO ijk = 1, lot 1763 1764 a10(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) + & 1765 & qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i))) 1766 a20(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) - & 1767 & qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i))) 1768 b10(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) + & 1769 & qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i))) 1770 b20(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) - & 1771 & qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i))) 1772 a11(ijk) = sin72*(b(ib+i)+b(ie+i)) + sin36*(b(ic+i)+b(id+i)) 1773 a21(ijk) = sin36*(b(ib+i)+b(ie+i)) - sin72*(b(ic+i)+b(id+i)) 1774 b11(ijk) = sin72*(a(ib+i)-a(ie+i)) + sin36*(a(ic+i)-a(id+i)) 1775 b21(ijk) = sin36*(a(ib+i)-a(ie+i)) - sin72*(a(ic+i)-a(id+i)) 1776 1777 c(ja+j) = a(ia+i) + ((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i))) 1778 d(ja+j) = b(ia+i) + ((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i))) 1779 c(jb+j) = c1*(a10(ijk)-a11(ijk)) - s1*(b10(ijk)+b11(ijk)) 1780 d(jb+j) = s1*(a10(ijk)-a11(ijk)) + c1*(b10(ijk)+b11(ijk)) 1781 c(je+j) = c4*(a10(ijk)+a11(ijk)) - s4*(b10(ijk)-b11(ijk)) 1782 d(je+j) = s4*(a10(ijk)+a11(ijk)) + c4*(b10(ijk)-b11(ijk)) 1783 c(jc+j) = c2*(a20(ijk)-a21(ijk)) - s2*(b20(ijk)+b21(ijk)) 1784 d(jc+j) = s2*(a20(ijk)-a21(ijk)) + c2*(b20(ijk)+b21(ijk)) 1785 c(jd+j) = c3*(a20(ijk)+a21(ijk)) - s3*(b20(ijk)-b21(ijk)) 1786 d(jd+j) = s3*(a20(ijk)+a21(ijk)) + c3*(b20(ijk)-b21(ijk)) 1787 1788 i = i + inc3 1789 j = j + inc4 1790 END DO 1791 ibase = ibase + inc1 1792 jbase = jbase + inc2 1793 END DO 1794 ia = ia + iink 1795 ib = ib + iink 1796 ic = ic + iink 1797 id = id - iink 1798 ie = ie - iink 1799 jbase = jbase + jump 1800 END DO 1801 IF (ib>id) GO TO 170 1802 110 CONTINUE 1803 ibase = 0 1804 DO l = 1, la 1805 i = ibase 1806 j = jbase 1807 !DIR$ IVDEP 1808 !CDIR NODEP 1809 !OCL NOVREC 1810 DO ijk = 1, lot 1811 c(ja+j) = (a(ia+i)+a(ib+i)) + a(ic+i) 1812 c(jb+j) = (qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - & 1813 & (sin36*b(ia+i)+sin72*b(ib+i)) 1814 c(je+j) = -(qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - & 1815 & (sin36*b(ia+i)+sin72*b(ib+i)) 1816 c(jc+j) = (qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - & 1817 & (sin72*b(ia+i)-sin36*b(ib+i)) 1818 c(jd+j) = -(qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - & 1819 & (sin72*b(ia+i)-sin36*b(ib+i)) 1820 i = i + inc3 1821 j = j + inc4 1822 END DO 1823 ibase = ibase + inc1 1824 jbase = jbase + inc2 1825 END DO 1826 1827 GO TO 170 1828 120 CONTINUE 1829 qqrt5 = 2.0_wp*qrt5 1830 ssin36 = 2.0_wp*sin36 1831 ssin72 = 2.0_wp*sin72 1832 DO l = 1, la 1833 i = ibase 1834 j = jbase 1835 !DIR$ IVDEP 1836 !CDIR NODEP 1837 !OCL NOVREC 1838 DO ijk = 1, lot 1839 c(ja+j) = 2.0_wp*(a(ia+i)+(a(ib+i)+a(ic+i))) 1840 c(jb+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ & 1841 & i))) - (ssin72*b(ib+i)+ssin36*b(ic+i)) 1842 c(jc+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ & 1843 & i))) - (ssin36*b(ib+i)-ssin72*b(ic+i)) 1844 c(jd+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ & 1845 & i))) + (ssin36*b(ib+i)-ssin72*b(ic+i)) 1846 c(je+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ & 1847 & i))) + (ssin72*b(ib+i)+ssin36*b(ic+i)) 1848 i = i + inc3 1849 j = j + inc4 1850 END DO 1851 ibase = ibase + inc1 1852 jbase = jbase + inc2 1853 END DO 1854 GO TO 170 1855 1856 ! Coding for factor 6 1857 1858 130 CONTINUE 1859 ia = 1 1860 ib = ia + (2*m-la)*inc1 1861 ic = ib + 2*m*inc1 1862 id = ic + 2*m*inc1 1863 ie = ic 1864 if = ib 1865 ja = 1 1866 jb = ja + jink 1867 jc = jb + jink 1868 jd = jc + jink 1869 je = jd + jink 1870 jf = je + jink 1871 1872 IF (la==m) GO TO 150 1873 1874 DO l = 1, la 1875 i = ibase 1876 j = jbase 1877 !DIR$ IVDEP 1878 !CDIR NODEP 1879 !OCL NOVREC 1880 DO ijk = 1, lot 1881 c(ja+j) = (a(ia+i)+a(id+i)) + (a(ib+i)+a(ic+i)) 1882 c(jd+j) = (a(ia+i)-a(id+i)) - (a(ib+i)-a(ic+i)) 1883 c(jb+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) - (sin60*(b(ib+ & 1884 & i)+b(ic+i))) 1885 c(jf+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) + (sin60*(b(ib+ & 1886 & i)+b(ic+i))) 1887 c(jc+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) - (sin60*(b(ib+ & 1888 & i)-b(ic+i))) 1889 c(je+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) + (sin60*(b(ib+ & 1890 & i)-b(ic+i))) 1891 i = i + inc3 1892 j = j + inc4 1893 END DO 1894 ibase = ibase + inc1 1895 jbase = jbase + inc2 1896 END DO 1897 ia = ia + iink 1898 iink = 2*iink 1899 ib = ib + iink 1900 ic = ic + iink 1901 id = id - iink 1902 ie = ie - iink 1903 if = if - iink 1904 jbase = jbase + jump 1905 jump = 2*jump + jink 1906 IF (ic==id) GO TO 140 1907 DO k = la, kstop, la 1908 kb = k + k 1909 kc = kb + kb 1910 kd = kc + kb 1911 ke = kd + kb 1912 kf = ke + kb 1913 c1 = trigs(kb+1) 1914 s1 = trigs(kb+2) 1915 c2 = trigs(kc+1) 1916 s2 = trigs(kc+2) 1917 c3 = trigs(kd+1) 1918 s3 = trigs(kd+2) 1919 c4 = trigs(ke+1) 1920 s4 = trigs(ke+2) 1921 c5 = trigs(kf+1) 1922 s5 = trigs(kf+2) 1923 ibase = 0 1924 DO l = 1, la 1925 i = ibase 1926 j = jbase 1927 !DIR$ IVDEP 1928 !CDIR NODEP 1929 !OCL NOVREC 1930 DO ijk = 1, lot 1931 1932 a11(ijk) = (a(ie+i)+a(ib+i)) + (a(ic+i)+a(if+i)) 1933 a20(ijk) = (a(ia+i)+a(id+i)) - 0.5_wp*a11(ijk) 1934 a21(ijk) = sin60*((a(ie+i)+a(ib+i))-(a(ic+i)+a(if+i))) 1935 b11(ijk) = (b(ib+i)-b(ie+i)) + (b(ic+i)-b(if+i)) 1936 b20(ijk) = (b(ia+i)-b(id+i)) - 0.5_wp*b11(ijk) 1937 b21(ijk) = sin60*((b(ib+i)-b(ie+i))-(b(ic+i)-b(if+i))) 1938 1939 c(ja+j) = (a(ia+i)+a(id+i)) + a11(ijk) 1940 d(ja+j) = (b(ia+i)-b(id+i)) + b11(ijk) 1941 c(jc+j) = c2*(a20(ijk)-b21(ijk)) - s2*(b20(ijk)+a21(ijk)) 1942 d(jc+j) = s2*(a20(ijk)-b21(ijk)) + c2*(b20(ijk)+a21(ijk)) 1943 c(je+j) = c4*(a20(ijk)+b21(ijk)) - s4*(b20(ijk)-a21(ijk)) 1944 d(je+j) = s4*(a20(ijk)+b21(ijk)) + c4*(b20(ijk)-a21(ijk)) 1945 1946 a11(ijk) = (a(ie+i)-a(ib+i)) + (a(ic+i)-a(if+i)) 1947 b11(ijk) = (b(ie+i)+b(ib+i)) - (b(ic+i)+b(if+i)) 1948 a20(ijk) = (a(ia+i)-a(id+i)) - 0.5_wp*a11(ijk) 1949 a21(ijk) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i))) 1950 b20(ijk) = (b(ia+i)+b(id+i)) + 0.5_wp*b11(ijk) 1951 b21(ijk) = sin60*((b(ie+i)+b(ib+i))+(b(ic+i)+b(if+i))) 1952 1953 c(jd+j) = c3*((a(ia+i)-a(id+i))+a11(ijk)) - s3*((b(ia+i)+b(id+ & 1954 & i))-b11(ijk)) 1955 d(jd+j) = s3*((a(ia+i)-a(id+i))+a11(ijk)) + c3*((b(ia+i)+b(id+ & 1956 & i))-b11(ijk)) 1957 c(jb+j) = c1*(a20(ijk)-b21(ijk)) - s1*(b20(ijk)-a21(ijk)) 1958 d(jb+j) = s1*(a20(ijk)-b21(ijk)) + c1*(b20(ijk)-a21(ijk)) 1959 c(jf+j) = c5*(a20(ijk)+b21(ijk)) - s5*(b20(ijk)+a21(ijk)) 1960 d(jf+j) = s5*(a20(ijk)+b21(ijk)) + c5*(b20(ijk)+a21(ijk)) 1961 1962 i = i + inc3 1963 j = j + inc4 1964 END DO 1965 ibase = ibase + inc1 1966 jbase = jbase + inc2 1967 END DO 1968 ia = ia + iink 1969 ib = ib + iink 1970 ic = ic + iink 1971 id = id - iink 1972 ie = ie - iink 1973 if = if - iink 1974 jbase = jbase + jump 1975 END DO 1976 IF (ic>id) GO TO 170 1977 140 CONTINUE 1978 ibase = 0 1979 DO l = 1, la 1980 i = ibase 1981 j = jbase 1982 !DIR$ IVDEP 1983 !CDIR NODEP 1984 !OCL NOVREC 1985 DO ijk = 1, lot 1986 c(ja+j) = a(ib+i) + (a(ia+i)+a(ic+i)) 1987 c(jd+j) = b(ib+i) - (b(ia+i)+b(ic+i)) 1988 c(jb+j) = (sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i)) 1989 c(jf+j) = -(sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i)) 1990 c(jc+j) = sin60*(b(ic+i)-b(ia+i)) + (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i)) 1991 c(je+j) = sin60*(b(ic+i)-b(ia+i)) - (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i)) 1992 i = i + inc3 1993 j = j + inc4 1994 END DO 1995 ibase = ibase + inc1 1996 jbase = jbase + inc2 1997 END DO 1998 1999 GO TO 170 2000 150 CONTINUE 2001 ssin60 = 2.0_wp*sin60 2002 DO l = 1, la 2003 i = ibase 2004 j = jbase 2005 !DIR$ IVDEP 2006 !CDIR NODEP 2007 !OCL NOVREC 2008 DO ijk = 1, lot 2009 c(ja+j) = (2.0_wp*(a(ia+i)+a(id+i))) + (2.0_wp*(a(ib+i)+a(ic+i))) 2010 c(jd+j) = (2.0_wp*(a(ia+i)-a(id+i))) - (2.0_wp*(a(ib+i)-a(ic+i))) 2011 c(jb+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) - (ssin60*(b(ib+ & 2012 & i)+b(ic+i))) 2013 c(jf+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) + (ssin60*(b(ib+ & 2014 & i)+b(ic+i))) 2015 c(jc+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) - (ssin60*(b(ib+ & 2016 & i)-b(ic+i))) 2017 c(je+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) + (ssin60*(b(ib+ & 2018 & i)-b(ic+i))) 2019 i = i + inc3 2020 j = j + inc4 2021 END DO 2022 ibase = ibase + inc1 2023 jbase = jbase + inc2 2024 END DO 2025 GO TO 170 2026 2027 ! Coding for factor 8 2028 2029 160 CONTINUE 2030 ibad = 3 2031 IF (la/=m) GO TO 180 2032 ia = 1 2033 ib = ia + la*inc1 2034 ic = ib + 2*la*inc1 2035 id = ic + 2*la*inc1 2036 ie = id + 2*la*inc1 2037 ja = 1 2038 jb = ja + jink 2039 jc = jb + jink 2040 jd = jc + jink 2041 je = jd + jink 2042 jf = je + jink 2043 jg = jf + jink 2044 jh = jg + jink 2045 ssin45 = SQRT(2.0_wp) 2046 2047 DO l = 1, la 2048 i = ibase 2049 j = jbase 2050 !DIR$ IVDEP 2051 !CDIR NODEP 2052 !OCL NOVREC 2053 DO ijk = 1, lot 2054 c(ja+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i))) 2055 c(je+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i))) 2056 c(jc+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i))) 2057 c(jg+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i))) 2058 c(jb+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) + ssin45*((a(ib+i)-a(id+ & 2059 & i))-(b(ib+i)+b(id+i))) 2060 c(jf+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) - ssin45*((a(ib+i)-a(id+ & 2061 & i))-(b(ib+i)+b(id+i))) 2062 c(jd+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) - ssin45*((a(ib+i)-a(id+ & 2063 & i))+(b(ib+i)+b(id+i))) 2064 c(jh+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) + ssin45*((a(ib+i)-a(id+ & 2065 & i))+(b(ib+i)+b(id+i))) 2066 i = i + inc3 2067 j = j + inc4 2068 END DO 2069 ibase = ibase + inc1 2070 jbase = jbase + inc2 2071 END DO 2072 2073 ! Return 2074 2075 170 CONTINUE 2076 ibad = 0 2077 180 CONTINUE 2078 ierr = ibad 2079 RETURN 2080 END SUBROUTINE rpassm 1387 IF ( igo == 7 ) igo = 6 1388 IF ( igo < 1 .OR. igo > 6 ) THEN 1389 ierr = 2 1390 RETURN 1391 ENDIF 1392 1393 1394 SELECT CASE ( igo ) 1395 ! 1396 !-- Coding for factor 2 1397 CASE ( 1 ) 1398 1399 ia = 1 1400 ib = ia + (2*m-la) * inc1 1401 ja = 1 1402 jb = ja + jink 1403 1404 IF ( la /= m ) THEN 1405 1406 DO l = 1, la 1407 i = ibase 1408 j = jbase 1409 !DIR$ IVDEP 1410 !CDIR NODEP 1411 !OCL NOVREC 1412 DO ijk = 1, lot 1413 c(ja+j) = a(ia+i) + a(ib+i) 1414 c(jb+j) = a(ia+i) - a(ib+i) 1415 i = i + inc3 1416 j = j + inc4 1417 ENDDO 1418 ibase = ibase + inc1 1419 jbase = jbase + inc2 1420 ENDDO 1421 1422 ia = ia + iink 1423 iink = 2 * iink 1424 ib = ib - iink 1425 ibase = 0 1426 jbase = jbase + jump 1427 jump = 2 * jump + jink 1428 1429 IF ( ia /= ib ) THEN 1430 1431 DO k = la, kstop, la 1432 kb = k + k 1433 c1 = trigs(kb+1) 1434 s1 = trigs(kb+2) 1435 ibase = 0 1436 DO l = 1, la 1437 i = ibase 1438 j = jbase 1439 !DIR$ IVDEP 1440 !CDIR NODEP 1441 !OCL NOVREC 1442 DO ijk = 1, lot 1443 c(ja+j) = a(ia+i) + a(ib+i) 1444 d(ja+j) = b(ia+i) - b(ib+i) 1445 c(jb+j) = c1*(a(ia+i)-a(ib+i)) - s1*(b(ia+i)+b(ib+i)) 1446 d(jb+j) = s1*(a(ia+i)-a(ib+i)) + c1*(b(ia+i)+b(ib+i)) 1447 i = i + inc3 1448 j = j + inc4 1449 ENDDO 1450 ibase = ibase + inc1 1451 jbase = jbase + inc2 1452 ENDDO 1453 ia = ia + iink 1454 ib = ib - iink 1455 jbase = jbase + jump 1456 ENDDO 1457 1458 IF ( ia > ib ) RETURN 1459 1460 ENDIF 1461 1462 ibase = 0 1463 DO l = 1, la 1464 i = ibase 1465 j = jbase 1466 !DIR$ IVDEP 1467 !CDIR NODEP 1468 !OCL NOVREC 1469 DO ijk = 1, lot 1470 c(ja+j) = a(ia+i) 1471 c(jb+j) = -b(ia+i) 1472 i = i + inc3 1473 j = j + inc4 1474 ENDDO 1475 ibase = ibase + inc1 1476 jbase = jbase + inc2 1477 ENDDO 1478 1479 ELSE 1480 1481 DO l = 1, la 1482 i = ibase 1483 j = jbase 1484 !DIR$ IVDEP 1485 !CDIR NODEP 1486 !OCL NOVREC 1487 DO ijk = 1, lot 1488 c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i)) 1489 c(jb+j) = 2.0_wp*(a(ia+i)-a(ib+i)) 1490 i = i + inc3 1491 j = j + inc4 1492 ENDDO 1493 ibase = ibase + inc1 1494 jbase = jbase + inc2 1495 ENDDO 1496 1497 ENDIF 1498 1499 ! 1500 !-- Coding for factor 3 1501 CASE ( 2 ) 1502 1503 ia = 1 1504 ib = ia + (2*m-la) * inc1 1505 ic = ib 1506 ja = 1 1507 jb = ja + jink 1508 jc = jb + jink 1509 1510 IF ( la /= m ) THEN 1511 1512 DO l = 1, la 1513 i = ibase 1514 j = jbase 1515 !DIR$ IVDEP 1516 !CDIR NODEP 1517 !OCL NOVREC 1518 DO ijk = 1, lot 1519 c(ja+j) = a(ia+i) + a(ib+i) 1520 c(jb+j) = (a(ia+i)-0.5_wp*a(ib+i)) - (sin60*(b(ib+i))) 1521 c(jc+j) = (a(ia+i)-0.5_wp*a(ib+i)) + (sin60*(b(ib+i))) 1522 i = i + inc3 1523 j = j + inc4 1524 ENDDO 1525 ibase = ibase + inc1 1526 jbase = jbase + inc2 1527 ENDDO 1528 ia = ia + iink 1529 iink = 2 * iink 1530 ib = ib + iink 1531 ic = ic - iink 1532 jbase = jbase + jump 1533 jump = 2 * jump + jink 1534 1535 IF ( ia /= ic ) THEN 1536 1537 DO k = la, kstop, la 1538 kb = k + k 1539 kc = kb + kb 1540 c1 = trigs(kb+1) 1541 s1 = trigs(kb+2) 1542 c2 = trigs(kc+1) 1543 s2 = trigs(kc+2) 1544 ibase = 0 1545 DO l = 1, la 1546 i = ibase 1547 j = jbase 1548 !DIR$ IVDEP 1549 !CDIR NODEP 1550 !OCL NOVREC 1551 DO ijk = 1, lot 1552 c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i)) 1553 d(ja+j) = b(ia+i) + (b(ib+i)-b(ic+i)) 1554 c(jb+j) = c1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ b(ic+i)))) & 1555 - s1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- a(ic+i)))) 1556 d(jb+j) = s1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ b(ic+i)))) & 1557 + c1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- a(ic+i)))) 1558 c(jc+j) = c2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ b(ic+i)))) & 1559 - s2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- a(ic+i)))) 1560 d(jc+j) = s2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ b(ic+i)))) & 1561 + c2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- a(ic+i)))) 1562 i = i + inc3 1563 j = j + inc4 1564 ENDDO 1565 ibase = ibase + inc1 1566 jbase = jbase + inc2 1567 ENDDO 1568 ia = ia + iink 1569 ib = ib + iink 1570 ic = ic - iink 1571 jbase = jbase + jump 1572 ENDDO 1573 1574 IF ( ia > ic ) RETURN 1575 1576 ENDIF 1577 1578 ibase = 0 1579 DO l = 1, la 1580 i = ibase 1581 j = jbase 1582 !DIR$ IVDEP 1583 !CDIR NODEP 1584 !OCL NOVREC 1585 DO ijk = 1, lot 1586 c(ja+j) = a(ia+i) + a(ib+i) 1587 c(jb+j) = (0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i)) 1588 c(jc+j) = -(0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i)) 1589 i = i + inc3 1590 j = j + inc4 1591 ENDDO 1592 ibase = ibase + inc1 1593 jbase = jbase + inc2 1594 ENDDO 1595 1596 ELSE 1597 1598 ssin60 = 2.0_wp * sin60 1599 DO l = 1, la 1600 i = ibase 1601 j = jbase 1602 !DIR$ IVDEP 1603 !CDIR NODEP 1604 !OCL NOVREC 1605 DO ijk = 1, lot 1606 c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i)) 1607 c(jb+j) = (2.0_wp*a(ia+i)-a(ib+i)) - (ssin60*b(ib+i)) 1608 c(jc+j) = (2.0_wp*a(ia+i)-a(ib+i)) + (ssin60*b(ib+i)) 1609 i = i + inc3 1610 j = j + inc4 1611 ENDDO 1612 ibase = ibase + inc1 1613 jbase = jbase + inc2 1614 ENDDO 1615 1616 ENDIF 1617 1618 ! 1619 !-- Coding for factor 4 1620 CASE ( 3 ) 1621 1622 ia = 1 1623 ib = ia + (2*m-la) * inc1 1624 ic = ib + 2*m*inc1 1625 id = ib 1626 ja = 1 1627 jb = ja + jink 1628 jc = jb + jink 1629 jd = jc + jink 1630 1631 IF ( la /= m ) THEN 1632 1633 DO l = 1, la 1634 i = ibase 1635 j = jbase 1636 !DIR$ IVDEP 1637 !CDIR NODEP 1638 !OCL NOVREC 1639 DO ijk = 1, lot 1640 c(ja+j) = (a(ia+i)+a(ic+i)) + a(ib+i) 1641 c(jb+j) = (a(ia+i)-a(ic+i)) - b(ib+i) 1642 c(jc+j) = (a(ia+i)+a(ic+i)) - a(ib+i) 1643 c(jd+j) = (a(ia+i)-a(ic+i)) + b(ib+i) 1644 i = i + inc3 1645 j = j + inc4 1646 ENDDO 1647 ibase = ibase + inc1 1648 jbase = jbase + inc2 1649 ENDDO 1650 ia = ia + iink 1651 iink = 2 * iink 1652 ib = ib + iink 1653 ic = ic - iink 1654 id = id - iink 1655 jbase = jbase + jump 1656 jump = 2 * jump + jink 1657 1658 IF ( ib /= ic ) THEN 1659 1660 DO k = la, kstop, la 1661 kb = k + k 1662 kc = kb + kb 1663 kd = kc + kb 1664 c1 = trigs(kb+1) 1665 s1 = trigs(kb+2) 1666 c2 = trigs(kc+1) 1667 s2 = trigs(kc+2) 1668 c3 = trigs(kd+1) 1669 s3 = trigs(kd+2) 1670 ibase = 0 1671 DO l = 1, la 1672 i = ibase 1673 j = jbase 1674 !DIR$ IVDEP 1675 !CDIR NODEP 1676 !OCL NOVREC 1677 DO ijk = 1, lot 1678 c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i)) 1679 d(ja+j) = (b(ia+i)-b(ic+i)) + (b(ib+i)-b(id+i)) 1680 c(jc+j) = c2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) - s2*((b(ia+i)-b(ic+i))& 1681 -(b(ib+i)-b(id+i))) 1682 d(jc+j) = s2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) + c2*((b(ia+i)-b(ic+i))& 1683 -(b(ib+i)-b(id+i))) 1684 c(jb+j) = c1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) - s1*((b(ia+i)+b(ic+i))& 1685 +(a(ib+i)-a(id+i))) 1686 d(jb+j) = s1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) + c1*((b(ia+i)+b(ic+i))& 1687 +(a(ib+i)-a(id+i))) 1688 c(jd+j) = c3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) - s3*((b(ia+i)+b(ic+i))& 1689 -(a(ib+i)-a(id+i))) 1690 d(jd+j) = s3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) + c3*((b(ia+i)+b(ic+i))& 1691 -(a(ib+i)-a(id+i))) 1692 i = i + inc3 1693 j = j + inc4 1694 ENDDO 1695 ibase = ibase + inc1 1696 jbase = jbase + inc2 1697 ENDDO 1698 ia = ia + iink 1699 ib = ib + iink 1700 ic = ic - iink 1701 id = id - iink 1702 jbase = jbase + jump 1703 ENDDO 1704 1705 IF ( ib > ic ) RETURN 1706 1707 ENDIF 1708 1709 ibase = 0 1710 sin45 = SQRT( 0.5_wp ) 1711 DO l = 1, la 1712 i = ibase 1713 j = jbase 1714 !DIR$ IVDEP 1715 !CDIR NODEP 1716 !OCL NOVREC 1717 DO ijk = 1, lot 1718 c(ja+j) = a(ia+i) + a(ib+i) 1719 c(jb+j) = sin45*((a(ia+i)-a(ib+i))-(b(ia+i)+b(ib+i))) 1720 c(jc+j) = b(ib+i) - b(ia+i) 1721 c(jd+j) = -sin45*((a(ia+i)-a(ib+i))+(b(ia+i)+b(ib+i))) 1722 i = i + inc3 1723 j = j + inc4 1724 ENDDO 1725 ibase = ibase + inc1 1726 jbase = jbase + inc2 1727 ENDDO 1728 1729 ELSE 1730 1731 DO l = 1, la 1732 i = ibase 1733 j = jbase 1734 !DIR$ IVDEP 1735 !CDIR NODEP 1736 !OCL NOVREC 1737 DO ijk = 1, lot 1738 c(ja+j) = 2.0_wp*((a(ia+i)+a(ic+i))+a(ib+i)) 1739 c(jb+j) = 2.0_wp*((a(ia+i)-a(ic+i))-b(ib+i)) 1740 c(jc+j) = 2.0_wp*((a(ia+i)+a(ic+i))-a(ib+i)) 1741 c(jd+j) = 2.0_wp*((a(ia+i)-a(ic+i))+b(ib+i)) 1742 i = i + inc3 1743 j = j + inc4 1744 ENDDO 1745 ibase = ibase + inc1 1746 jbase = jbase + inc2 1747 ENDDO 1748 1749 ENDIF 1750 1751 ! 1752 !-- Coding for factor 5 1753 CASE ( 4 ) 1754 1755 ia = 1 1756 ib = ia + (2*m-la) * inc1 1757 ic = ib + 2*m*inc1 1758 id = ic 1759 ie = ib 1760 ja = 1 1761 jb = ja + jink 1762 jc = jb + jink 1763 jd = jc + jink 1764 je = jd + jink 1765 1766 IF ( la /= m ) THEN 1767 1768 DO l = 1, la 1769 i = ibase 1770 j = jbase 1771 !DIR$ IVDEP 1772 !CDIR NODEP 1773 !OCL NOVREC 1774 DO ijk = 1, lot 1775 c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i)) 1776 c(jb+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) - & 1777 (sin72*b(ib+i)+sin36*b(ic+i)) 1778 c(jc+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) - & 1779 (sin36*b(ib+i)-sin72*b(ic+i)) 1780 c(jd+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) + & 1781 (sin36*b(ib+i)-sin72*b(ic+i)) 1782 c(je+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) + & 1783 (sin72*b(ib+i)+sin36*b(ic+i)) 1784 i = i + inc3 1785 j = j + inc4 1786 ENDDO 1787 ibase = ibase + inc1 1788 jbase = jbase + inc2 1789 ENDDO 1790 ia = ia + iink 1791 iink = 2 * iink 1792 ib = ib + iink 1793 ic = ic + iink 1794 id = id - iink 1795 ie = ie - iink 1796 jbase = jbase + jump 1797 jump = 2 * jump + jink 1798 1799 IF ( ib /= id ) THEN 1800 1801 DO k = la, kstop, la 1802 kb = k + k 1803 kc = kb + kb 1804 kd = kc + kb 1805 ke = kd + kb 1806 c1 = trigs(kb+1) 1807 s1 = trigs(kb+2) 1808 c2 = trigs(kc+1) 1809 s2 = trigs(kc+2) 1810 c3 = trigs(kd+1) 1811 s3 = trigs(kd+2) 1812 c4 = trigs(ke+1) 1813 s4 = trigs(ke+2) 1814 ibase = 0 1815 DO l = 1, la 1816 i = ibase 1817 j = jbase 1818 !DIR$ IVDEP 1819 !CDIR NODEP 1820 !OCL NOVREC 1821 DO ijk = 1, lot 1822 a10(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) + & 1823 qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i))) 1824 a20(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) - & 1825 qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i))) 1826 b10(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) + & 1827 qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i))) 1828 b20(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) - & 1829 qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i))) 1830 a11(ijk) = sin72*(b(ib+i)+b(ie+i)) + sin36*(b(ic+i)+b(id+i)) 1831 a21(ijk) = sin36*(b(ib+i)+b(ie+i)) - sin72*(b(ic+i)+b(id+i)) 1832 b11(ijk) = sin72*(a(ib+i)-a(ie+i)) + sin36*(a(ic+i)-a(id+i)) 1833 b21(ijk) = sin36*(a(ib+i)-a(ie+i)) - sin72*(a(ic+i)-a(id+i)) 1834 c(ja+j) = a(ia+i) + ((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i))) 1835 d(ja+j) = b(ia+i) + ((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i))) 1836 c(jb+j) = c1*(a10(ijk)-a11(ijk)) - s1*(b10(ijk)+b11(ijk)) 1837 d(jb+j) = s1*(a10(ijk)-a11(ijk)) + c1*(b10(ijk)+b11(ijk)) 1838 c(je+j) = c4*(a10(ijk)+a11(ijk)) - s4*(b10(ijk)-b11(ijk)) 1839 d(je+j) = s4*(a10(ijk)+a11(ijk)) + c4*(b10(ijk)-b11(ijk)) 1840 c(jc+j) = c2*(a20(ijk)-a21(ijk)) - s2*(b20(ijk)+b21(ijk)) 1841 d(jc+j) = s2*(a20(ijk)-a21(ijk)) + c2*(b20(ijk)+b21(ijk)) 1842 c(jd+j) = c3*(a20(ijk)+a21(ijk)) - s3*(b20(ijk)-b21(ijk)) 1843 d(jd+j) = s3*(a20(ijk)+a21(ijk)) + c3*(b20(ijk)-b21(ijk)) 1844 1845 i = i + inc3 1846 j = j + inc4 1847 ENDDO 1848 ibase = ibase + inc1 1849 jbase = jbase + inc2 1850 ENDDO 1851 ia = ia + iink 1852 ib = ib + iink 1853 ic = ic + iink 1854 id = id - iink 1855 ie = ie - iink 1856 jbase = jbase + jump 1857 ENDDO 1858 1859 IF ( ib > id ) RETURN 1860 1861 ENDIF 1862 1863 ibase = 0 1864 DO l = 1, la 1865 i = ibase 1866 j = jbase 1867 !DIR$ IVDEP 1868 !CDIR NODEP 1869 !OCL NOVREC 1870 DO ijk = 1, lot 1871 c(ja+j) = (a(ia+i)+a(ib+i)) + a(ic+i) 1872 c(jb+j) = (qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - & 1873 (sin36*b(ia+i)+sin72*b(ib+i)) 1874 c(je+j) = -(qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - & 1875 (sin36*b(ia+i)+sin72*b(ib+i)) 1876 c(jc+j) = (qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - & 1877 (sin72*b(ia+i)-sin36*b(ib+i)) 1878 c(jd+j) = -(qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - & 1879 (sin72*b(ia+i)-sin36*b(ib+i)) 1880 i = i + inc3 1881 j = j + inc4 1882 ENDDO 1883 ibase = ibase + inc1 1884 jbase = jbase + inc2 1885 ENDDO 1886 1887 ELSE 1888 1889 qqrt5 = 2.0_wp * qrt5 1890 ssin36 = 2.0_wp * sin36 1891 ssin72 = 2.0_wp * sin72 1892 DO l = 1, la 1893 i = ibase 1894 j = jbase 1895 !DIR$ IVDEP 1896 !CDIR NODEP 1897 !OCL NOVREC 1898 DO ijk = 1, lot 1899 c(ja+j) = 2.0_wp*(a(ia+i)+(a(ib+i)+a(ic+i))) 1900 c(jb+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+i))) & 1901 - (ssin72*b(ib+i)+ssin36*b(ic+i)) 1902 c(jc+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+i))) & 1903 - (ssin36*b(ib+i)-ssin72*b(ic+i)) 1904 c(jd+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+i))) & 1905 + (ssin36*b(ib+i)-ssin72*b(ic+i)) 1906 c(je+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+i))) & 1907 + (ssin72*b(ib+i)+ssin36*b(ic+i)) 1908 i = i + inc3 1909 j = j + inc4 1910 ENDDO 1911 ibase = ibase + inc1 1912 jbase = jbase + inc2 1913 ENDDO 1914 1915 ENDIF 1916 1917 ! 1918 !-- Coding for factor 6 1919 CASE ( 5 ) 1920 1921 ia = 1 1922 ib = ia + (2*m-la) * inc1 1923 ic = ib + 2*m*inc1 1924 id = ic + 2*m*inc1 1925 ie = ic 1926 if = ib 1927 ja = 1 1928 jb = ja + jink 1929 jc = jb + jink 1930 jd = jc + jink 1931 je = jd + jink 1932 jf = je + jink 1933 1934 IF ( la /= m ) THEN 1935 1936 DO l = 1, la 1937 i = ibase 1938 j = jbase 1939 !DIR$ IVDEP 1940 !CDIR NODEP 1941 !OCL NOVREC 1942 DO ijk = 1, lot 1943 c(ja+j) = (a(ia+i)+a(id+i)) + (a(ib+i)+a(ic+i)) 1944 c(jd+j) = (a(ia+i)-a(id+i)) - (a(ib+i)-a(ic+i)) 1945 c(jb+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) - (sin60*(b(ib+i)+b(ic+i))) 1946 c(jf+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) + (sin60*(b(ib+i)+b(ic+i))) 1947 c(jc+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) - (sin60*(b(ib+i)-b(ic+i))) 1948 c(je+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) + (sin60*(b(ib+i)-b(ic+i))) 1949 i = i + inc3 1950 j = j + inc4 1951 ENDDO 1952 ibase = ibase + inc1 1953 jbase = jbase + inc2 1954 ENDDO 1955 ia = ia + iink 1956 iink = 2 * iink 1957 ib = ib + iink 1958 ic = ic + iink 1959 id = id - iink 1960 ie = ie - iink 1961 if = if - iink 1962 jbase = jbase + jump 1963 jump = 2 * jump + jink 1964 1965 IF ( ic /= id ) THEN 1966 1967 DO k = la, kstop, la 1968 kb = k + k 1969 kc = kb + kb 1970 kd = kc + kb 1971 ke = kd + kb 1972 kf = ke + kb 1973 c1 = trigs(kb+1) 1974 s1 = trigs(kb+2) 1975 c2 = trigs(kc+1) 1976 s2 = trigs(kc+2) 1977 c3 = trigs(kd+1) 1978 s3 = trigs(kd+2) 1979 c4 = trigs(ke+1) 1980 s4 = trigs(ke+2) 1981 c5 = trigs(kf+1) 1982 s5 = trigs(kf+2) 1983 ibase = 0 1984 DO l = 1, la 1985 i = ibase 1986 j = jbase 1987 !DIR$ IVDEP 1988 !CDIR NODEP 1989 !OCL NOVREC 1990 DO ijk = 1, lot 1991 a11(ijk) = (a(ie+i)+a(ib+i)) + (a(ic+i)+a(if+i)) 1992 a20(ijk) = (a(ia+i)+a(id+i)) - 0.5_wp*a11(ijk) 1993 a21(ijk) = sin60*((a(ie+i)+a(ib+i))-(a(ic+i)+a(if+i))) 1994 b11(ijk) = (b(ib+i)-b(ie+i)) + (b(ic+i)-b(if+i)) 1995 b20(ijk) = (b(ia+i)-b(id+i)) - 0.5_wp*b11(ijk) 1996 b21(ijk) = sin60*((b(ib+i)-b(ie+i))-(b(ic+i)-b(if+i))) 1997 1998 c(ja+j) = (a(ia+i)+a(id+i)) + a11(ijk) 1999 d(ja+j) = (b(ia+i)-b(id+i)) + b11(ijk) 2000 c(jc+j) = c2*(a20(ijk)-b21(ijk)) - s2*(b20(ijk)+a21(ijk)) 2001 d(jc+j) = s2*(a20(ijk)-b21(ijk)) + c2*(b20(ijk)+a21(ijk)) 2002 c(je+j) = c4*(a20(ijk)+b21(ijk)) - s4*(b20(ijk)-a21(ijk)) 2003 d(je+j) = s4*(a20(ijk)+b21(ijk)) + c4*(b20(ijk)-a21(ijk)) 2004 2005 a11(ijk) = (a(ie+i)-a(ib+i)) + (a(ic+i)-a(if+i)) 2006 b11(ijk) = (b(ie+i)+b(ib+i)) - (b(ic+i)+b(if+i)) 2007 a20(ijk) = (a(ia+i)-a(id+i)) - 0.5_wp*a11(ijk) 2008 a21(ijk) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i))) 2009 b20(ijk) = (b(ia+i)+b(id+i)) + 0.5_wp*b11(ijk) 2010 b21(ijk) = sin60*((b(ie+i)+b(ib+i))+(b(ic+i)+b(if+i))) 2011 2012 c(jd+j) = c3*((a(ia+i)-a(id+i))+a11(ijk)) - s3*((b(ia+i)+b(id+i))-b11(ijk)) 2013 d(jd+j) = s3*((a(ia+i)-a(id+i))+a11(ijk)) + c3*((b(ia+i)+b(id+i))-b11(ijk)) 2014 c(jb+j) = c1*(a20(ijk)-b21(ijk)) - s1*(b20(ijk)-a21(ijk)) 2015 d(jb+j) = s1*(a20(ijk)-b21(ijk)) + c1*(b20(ijk)-a21(ijk)) 2016 c(jf+j) = c5*(a20(ijk)+b21(ijk)) - s5*(b20(ijk)+a21(ijk)) 2017 d(jf+j) = s5*(a20(ijk)+b21(ijk)) + c5*(b20(ijk)+a21(ijk)) 2018 2019 i = i + inc3 2020 j = j + inc4 2021 ENDDO 2022 ibase = ibase + inc1 2023 jbase = jbase + inc2 2024 ENDDO 2025 ia = ia + iink 2026 ib = ib + iink 2027 ic = ic + iink 2028 id = id - iink 2029 ie = ie - iink 2030 if = if - iink 2031 jbase = jbase + jump 2032 ENDDO 2033 2034 IF ( ic > id ) RETURN 2035 2036 ENDIF 2037 2038 ibase = 0 2039 DO l = 1, la 2040 i = ibase 2041 j = jbase 2042 !DIR$ IVDEP 2043 !CDIR NODEP 2044 !OCL NOVREC 2045 DO ijk = 1, lot 2046 c(ja+j) = a(ib+i) + (a(ia+i)+a(ic+i)) 2047 c(jd+j) = b(ib+i) - (b(ia+i)+b(ic+i)) 2048 c(jb+j) = (sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i)) 2049 c(jf+j) = -(sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i)) 2050 c(jc+j) = sin60*(b(ic+i)-b(ia+i)) + (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i)) 2051 c(je+j) = sin60*(b(ic+i)-b(ia+i)) - (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i)) 2052 i = i + inc3 2053 j = j + inc4 2054 ENDDO 2055 ibase = ibase + inc1 2056 jbase = jbase + inc2 2057 ENDDO 2058 2059 ELSE 2060 2061 ssin60 = 2.0_wp * sin60 2062 DO l = 1, la 2063 i = ibase 2064 j = jbase 2065 !DIR$ IVDEP 2066 !CDIR NODEP 2067 !OCL NOVREC 2068 DO ijk = 1, lot 2069 c(ja+j) = (2.0_wp*(a(ia+i)+a(id+i))) + (2.0_wp*(a(ib+i)+a(ic+i))) 2070 c(jd+j) = (2.0_wp*(a(ia+i)-a(id+i))) - (2.0_wp*(a(ib+i)-a(ic+i))) 2071 c(jb+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) - (ssin60*(b(ib+i)+b(ic+i))) 2072 c(jf+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) + (ssin60*(b(ib+i)+b(ic+i))) 2073 c(jc+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) - (ssin60*(b(ib+i)-b(ic+i))) 2074 c(je+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) + (ssin60*(b(ib+i)-b(ic+i))) 2075 i = i + inc3 2076 j = j + inc4 2077 ENDDO 2078 ibase = ibase + inc1 2079 jbase = jbase + inc2 2080 ENDDO 2081 2082 ENDIF 2083 2084 ! 2085 !-- Coding for factor 8 2086 CASE ( 6 ) 2087 2088 IF ( la /= m ) THEN 2089 ierr = 3 2090 RETURN 2091 ENDIF 2092 2093 ia = 1 2094 ib = ia + la*inc1 2095 ic = ib + 2*la*inc1 2096 id = ic + 2*la*inc1 2097 ie = id + 2*la*inc1 2098 ja = 1 2099 jb = ja + jink 2100 jc = jb + jink 2101 jd = jc + jink 2102 je = jd + jink 2103 jf = je + jink 2104 jg = jf + jink 2105 jh = jg + jink 2106 ssin45 = SQRT( 2.0_wp ) 2107 2108 DO l = 1, la 2109 i = ibase 2110 j = jbase 2111 !DIR$ IVDEP 2112 !CDIR NODEP 2113 !OCL NOVREC 2114 DO ijk = 1, lot 2115 c(ja+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i))) 2116 c(je+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i))) 2117 c(jc+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i))) 2118 c(jg+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i))) 2119 c(jb+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) + ssin45*((a(ib+i)-a(id+i))-(b(ib+i)+b(id+i))) 2120 c(jf+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) - ssin45*((a(ib+i)-a(id+i))-(b(ib+i)+b(id+i))) 2121 c(jd+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) - ssin45*((a(ib+i)-a(id+i))+(b(ib+i)+b(id+i))) 2122 c(jh+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) + ssin45*((a(ib+i)-a(id+i))+(b(ib+i)+b(id+i))) 2123 i = i + inc3 2124 j = j + inc4 2125 ENDDO 2126 ibase = ibase + inc1 2127 jbase = jbase + inc2 2128 ENDDO 2129 2130 END SELECT 2131 2132 END SUBROUTINE rpassm 2081 2133 2082 2134 !------------------------------------------------------------------------------! … … 2088 2140 !> functins required by fft99 & fft991cy 2089 2141 !------------------------------------------------------------------------------! 2090 SUBROUTINE set99(trigs,ifax,n)2091 2092 2093 USE control_parameters, &2142 SUBROUTINE set99( trigs, ifax, n ) 2143 2144 2145 USE control_parameters, & 2094 2146 ONLY: message_string 2095 2147 … … 2098 2150 IMPLICIT NONE 2099 2151 2100 ! Scalar arguments2101 2152 INTEGER(iwp) :: n !< 2102 2153 2103 ! Array arguments2104 2154 INTEGER(iwp) :: ifax(*) !< 2105 2155 REAL(wp) :: trigs(*) !< 2106 2156 2107 2157 2108 ! Local scalars:2109 2158 REAL(wp) :: angle !< 2110 2159 REAL(wp) :: del !< … … 2119 2168 INTEGER(iwp) :: nu !< 2120 2169 2121 ! Local arrays:2122 2170 INTEGER(iwp) :: jfax(10) !< 2123 2171 INTEGER(iwp) :: lfax(7) !< 2124 2172 2125 ! Intrinsic functions 2126 ! INTRINSIC ASIN, COS, MOD, REAL, SIN 2127 2128 ! Data statements 2129 DATA lfax/6, 8, 5, 4, 3, 2, 1/ 2130 2131 2132 ! Executable statements 2173 DATA lfax / 6, 8, 5, 4, 3, 2, 1 / 2174 2175 2133 2176 ixxx = 1 2134 2177 2135 del = 4.0_wp *ASIN(1.0_wp)/REAL(n,KIND=wp)2178 del = 4.0_wp * ASIN( 1.0_wp ) / REAL( n, KIND=wp ) 2136 2179 nil = 0 2137 2180 nhl = (n/2) - 1 2138 DO k = nil, nhl 2139 angle = REAL(k,KIND=wp)*del 2140 trigs(2*k+1) = COS(angle) 2141 trigs(2*k+2) = SIN(angle) 2142 END DO 2143 2144 ! Find factors of n (8,6,5,4,3,2; only one 8 allowed) 2145 ! Look for sixes first, store factors in descending order 2146 nu = n 2181 DO k = nil, nhl 2182 angle = REAL( k, KIND=wp ) * del 2183 trigs(2*k+1) = COS( angle ) 2184 trigs(2*k+2) = SIN( angle ) 2185 ENDDO 2186 2187 ! 2188 !-- Find factors of n (8,6,5,4,3,2; only one 8 allowed). 2189 !-- Look for sixes first, store factors in descending order. 2190 nu = n 2147 2191 ifac = 6 2148 k = 0 2149 l = 1 2150 10 CONTINUE 2151 IF (MOD(nu,ifac)/=0) GO TO 30 2152 k = k + 1 2153 jfax(k) = ifac 2154 IF (ifac/=8) GO TO 20 2155 IF (k==1) GO TO 20 2156 jfax(1) = 8 2157 jfax(k) = 6 2158 20 CONTINUE 2159 nu = nu/ifac 2160 IF (nu==1) GO TO 40 2161 IF (ifac/=8) GO TO 10 2162 30 CONTINUE 2163 l = l + 1 2164 ifac = lfax(l) 2165 IF (ifac>1) GO TO 10 2166 2167 ! WRITE (nout,'(A,I4,A)') ' n =',n,' - Contains illegal factors' 2168 message_string = 'number of gridpoints along x or/and y ' // & 2169 'contain illegal factors' // & 2170 '&only factors 2, 3, 5 are allowed' 2171 CALL message( 'temperton_fft', 'PA0311', 1, 2, 0, 6, 0 ) 2172 2173 RETURN 2174 2175 ! Now reverse order of factors 2176 40 CONTINUE 2192 k = 0 2193 l = 1 2194 2195 DO 2196 2197 IF ( MOD( nu, ifac ) == 0 ) THEN 2198 k = k + 1 2199 jfax(k) = ifac 2200 IF ( ifac == 8 .AND. k /= 1 ) THEN 2201 jfax(1) = 8 2202 jfax(k) = 6 2203 ENDIF 2204 nu = nu/ifac 2205 IF ( nu == 1 ) EXIT 2206 IF ( ifac /= 8 ) CYCLE 2207 ENDIF 2208 2209 l = l + 1 2210 ifac = lfax(l) 2211 IF (ifac < 2 ) THEN 2212 message_string = 'number of gridpoints along x or/and y ' // & 2213 'contain illegal factors' // & 2214 '&only factors 2, 3, 5 are allowed' 2215 CALL message( 'temperton_fft', 'PA0311', 1, 2, 0, 6, 0 ) 2216 RETURN 2217 ENDIF 2218 2219 ENDDO 2220 ! 2221 !-- Now reverse order of factors 2177 2222 nfax = k 2178 2223 ifax(1) = nfax 2179 DO i = 1, nfax2224 DO i = 1, nfax 2180 2225 ifax(nfax+2-i) = jfax(i) 2181 END 2226 ENDDO 2182 2227 ifax(10) = n 2183 RETURN 2184 2228 2229 END SUBROUTINE set99 2185 2230 2186 2231 END MODULE temperton_fft -
palm/trunk/SOURCE/wind_turbine_model_mod.f90
r3685 r3725 26 26 ! ----------------- 27 27 ! $Id$ 28 ! unused variables removed 29 ! 30 ! 3685 2019-01-21 01:02:11Z knoop 28 31 ! Some interface calls moved to module_interface + cleanup 29 32 ! … … 237 240 REAL(wp), DIMENSION(1:100) :: rnac = 0.0_wp !< nacelle diameter [m] 238 241 REAL(wp), DIMENSION(1:100) :: rr = 63.0_wp !< rotor radius [m] 239 REAL(wp), DIMENSION(1:100) :: turb_cd_nacelle = 0.85_wp !< drag coefficient for nacelle242 ! REAL(wp), DIMENSION(1:100) :: turb_cd_nacelle = 0.85_wp !< drag coefficient for nacelle 240 243 REAL(wp), DIMENSION(1:100) :: turb_cd_tower = 1.2_wp !< drag coefficient for tower 241 244 … … 306 309 REAL(wp) :: eps_min !< 307 310 REAL(wp) :: eps_min2 !< 308 REAL(wp) :: sqrt_arg !<309 311 310 312 ! … … 519 521 IMPLICIT NONE 520 522 521 INTEGER(iwp) :: ierrn !<522 523 523 CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file 524 524 … … 531 531 rnac, rr, segment_length, segment_width, & 532 532 slope2, speed_control, tilt, time_turbine_on,& 533 turb_cd_ nacelle, turb_cd_tower, pitch_rate,&533 turb_cd_tower, pitch_rate, & 534 534 yaw_control, yaw_speed, tl_cor 535 ! , turb_cd_nacelle 535 536 536 537 NAMELIST /wind_turbine_parameters/ & … … 543 544 rnac, rr, segment_length, segment_width, & 544 545 slope2, speed_control, tilt, time_turbine_on,& 545 turb_cd_ nacelle, turb_cd_tower, pitch_rate,&546 turb_cd_tower, pitch_rate, & 546 547 yaw_control, yaw_speed, tl_cor 548 ! , turb_cd_nacelle 547 549 ! 548 550 !-- Try to find wind turbine model package … … 980 982 INTEGER(iwp) :: tower_n !< 981 983 INTEGER(iwp) :: tower_s !< 982 !983 !-- Help variables for the calulaction of the nacelle drag984 INTEGER(iwp) :: i_ip !<985 INTEGER(iwp) :: i_ipg !<986 987 REAL(wp) :: yvalue988 REAL(wp) :: dy_int !<989 REAL(wp) :: dz_int !<990 984 991 985 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: circle_points !<
Note: See TracChangeset
for help on using the changeset viewer.
