|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00016 MODULE pw_spline_utils 00017 USE cell_types, ONLY: cell_type 00018 USE cp_output_handling, ONLY: cp_add_iter_level,& 00019 cp_iterate,& 00020 cp_rm_iter_level 00021 USE input_constants, ONLY: no_precond,& 00022 precond_spl3_1,& 00023 precond_spl3_2,& 00024 precond_spl3_3,& 00025 precond_spl3_aint,& 00026 precond_spl3_aint2,& 00027 spline3_nopbc_interp,& 00028 spline3_pbc_interp 00029 USE input_section_types, ONLY: section_vals_type,& 00030 section_vals_val_get 00031 USE kinds, ONLY: dp 00032 USE mathconstants, ONLY: twopi 00033 USE message_passing, ONLY: mp_alltoall,& 00034 mp_comm_compare,& 00035 mp_sendrecv,& 00036 mp_sum 00037 USE pw_grid_types, ONLY: FULLSPACE,& 00038 PW_MODE_LOCAL 00039 USE pw_methods, ONLY: pw_axpy,& 00040 pw_copy,& 00041 pw_integral_aa_fast,& 00042 pw_integral_ab_fast,& 00043 pw_zero 00044 USE pw_pool_types, ONLY: pw_pool_create_pw,& 00045 pw_pool_give_back_pw,& 00046 pw_pool_release,& 00047 pw_pool_retain,& 00048 pw_pool_type 00049 USE pw_types, ONLY: COMPLEXDATA1D,& 00050 REALDATA3D,& 00051 REALSPACE,& 00052 RECIPROCALSPACE,& 00053 SQUARE,& 00054 pw_p_type,& 00055 pw_type 00056 USE timings, ONLY: timeset,& 00057 timestop 00058 #include "cp_common_uses.h" 00059 00060 IMPLICIT NONE 00061 PRIVATE 00062 00063 LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.TRUE. 00064 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_spline_utils' 00065 00066 REAL ( KIND = dp ), PUBLIC, PARAMETER, DIMENSION(4) :: nn10_coeffs= 00067 (/ 125._dp/216._dp, 25._dp/432._dp, 5._dp/864._dp, 1._dp/1728._dp /), 00068 spline3_coeffs= 00069 (/ 8._dp/(27._dp),2._dp/(27._dp),1._dp/(27._dp*2._dp), 00070 1._dp/(27._dp*8._dp) /), 00071 spline2_coeffs= 00072 (/ 27._dp/(64._dp),9._dp/(64._dp*2_dp),3._dp/(64._dp*4._dp), 00073 1._dp/(64._dp*8._dp) /), 00074 nn50_coeffs= 00075 (/ 15625._dp/17576._dp, 625._dp/35152._dp, 25._dp/70304._dp, 00076 1._dp/140608._dp /), 00077 spl3_aint_coeff= 00078 (/ 46._dp/27._dp,-2._dp/(27._dp),-1._dp/(27._dp*2._dp), 00079 -1._dp/(27._dp*8._dp) /), 00080 spl3_precond1_coeff= 00081 (/ 64._dp/3._dp,-8._dp/3._dp,-1._dp/3._dp, -1._dp/24._dp /), 00082 spl3_1d_transf_coeffs= 00083 (/ 2._dp/3._dp,23._dp/48._dp,1._dp/6._dp,1._dp/48._dp /) 00084 00085 REAL ( KIND = dp ), PUBLIC, PARAMETER, DIMENSION(3) :: spline3_deriv_coeffs= 00086 (/ 2.0_dp/9.0_dp, 1.0_dp/18.0_dp, 1.0_dp/72.0_dp /), 00087 spline2_deriv_coeffs= 00088 (/ 9.0_dp/32.0_dp, 3.0_dp/64.0_dp, 1.0_dp/128.0_dp /), 00089 nn10_deriv_coeffs= 00090 (/ 25._dp/72._dp, 5._dp/144, 1._dp/288._dp /), 00091 nn50_deriv_coeffs= 00092 (/ 625._dp/1352._dp, 25._dp/2704._dp,1._dp/5408._dp /), 00093 spl3_1d_coeffs0= 00094 (/ 1._dp/6_dp,2._dp/3._dp,1._dp/6._dp /), 00095 spl3_1d_transf_border1= 00096 (/ 0.517977704_dp, 0.464044595_dp, 0.17977701e-1_dp /) 00097 00098 INTEGER, SAVE, PRIVATE :: last_precond_id=0 00099 00100 PUBLIC :: pw_spline3_interpolate_values_g, pw_spline3_evaluate_values_g,& 00101 pw_spline3_deriv_g 00102 PUBLIC :: pw_spline_scale_deriv, pw_nn_smear_g 00103 PUBLIC :: pw_spline2_interpolate_values_g, pw_spline2_evaluate_values_g,& 00104 pw_spline2_deriv_g 00105 PUBLIC :: pw_nn_compose_r, pw_nn_smear_r, pw_nn_deriv_r, pw_prolongate_s3,& 00106 pw_restrict_s3, spl3_nopbc, spl3_pbc, pw_nn_compose_r_no_pbc, spl3_nopbct 00107 PUBLIC :: add_fine2coarse 00108 PUBLIC :: pw_spline_precond_create,& 00109 pw_spline_do_precond,& 00110 pw_spline_precond_set_kind,& 00111 find_coeffs,& 00112 pw_spline_precond_release,& 00113 pw_spline_precond_type,& 00114 Eval_Interp_Spl3_pbc,& 00115 Eval_d_Interp_Spl3_pbc 00116 !*** 00117 00118 ! ***************************************************************************** 00123 TYPE pw_spline_precond_type 00124 INTEGER :: ref_count,id_nr, kind 00125 REAL(kind=dp), DIMENSION(4) :: coeffs 00126 REAL(kind=dp), DIMENSION(3) :: coeffs_1d 00127 LOGICAL :: sharpen, normalize, pbc, transpose 00128 TYPE(pw_pool_type), POINTER :: pool 00129 END TYPE pw_spline_precond_type 00130 00131 CONTAINS 00132 00133 ! ***************************************************************************** 00146 SUBROUTINE pw_spline2_interpolate_values_g(spline_g,error) 00147 TYPE(pw_type), POINTER :: spline_g 00148 TYPE(cp_error_type), INTENT(inout) :: error 00149 00150 CHARACTER(len=*), PARAMETER :: 00151 routineN = 'pw_spline2_interpolate_values_g', 00152 routineP = moduleN//':'//routineN 00153 00154 INTEGER :: handle, i, ii, j, k, stat 00155 INTEGER, DIMENSION(2, 3) :: gbo 00156 INTEGER, DIMENSION(3) :: n_tot 00157 LOGICAL :: failure 00158 REAL(KIND=dp) :: c23, coeff 00159 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cosIVals, cosJVals, cosKVals 00160 00161 CALL timeset(routineN,handle) 00162 00163 failure=.FALSE. 00164 n_tot(1:3) = spline_g%pw_grid%npts (1:3) 00165 gbo = spline_g%pw_grid%bounds 00166 00167 CPPrecondition(ASSOCIATED(spline_g),cp_failure_level,routineP,error,failure) 00168 CPPrecondition(spline_g%in_use==COMPLEXDATA1D,cp_failure_level,routineP,error,failure) 00169 CPPrecondition(spline_g%in_space==RECIPROCALSPACE,cp_failure_level,routineP,error,failure) 00170 CPPrecondition(.NOT.spline_g%pw_grid%spherical,cp_failure_level,routineP,error,failure) 00171 CPPrecondition(spline_g%pw_grid%grid_span==FULLSPACE,cp_failure_level,routineP,error,failure) 00172 00173 ALLOCATE(cosIVals(gbo(1,1):gbo(2,1)),cosJVals(gbo(1,2):gbo(2,2)),& 00174 cosKVals(gbo(1,3):gbo(2,3)),stat=stat) 00175 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00176 00177 coeff=twopi/n_tot(1) 00178 !$omp parallel do private(i) 00179 DO i=gbo(1,1),gbo(2,1) 00180 cosIVals(i)=COS(coeff*REAL(i,dp)) 00181 END DO 00182 coeff=twopi/n_tot(2) 00183 !$omp parallel do private(j) 00184 DO j=gbo(1,2),gbo(2,2) 00185 cosJVals(j)=COS(coeff*REAL(j,dp)) 00186 END DO 00187 coeff=twopi/n_tot(3) 00188 !$omp parallel do private(k) 00189 DO k=gbo(1,3),gbo(2,3) 00190 cosKVals(k)=COS(coeff*REAL(k,dp)) 00191 END DO 00192 00193 !$omp parallel do private(i,j,k,ii,coeff,c23) 00194 DO ii=1,SIZE(spline_g%cc) 00195 i=spline_g%pw_grid%g_hat(1,ii) 00196 j=spline_g%pw_grid%g_hat(2,ii) 00197 k=spline_g%pw_grid%g_hat(3,ii) 00198 00199 c23=cosJVals(j)*cosKVals(k) 00200 coeff=64.0_dp/(cosIVals(i)*c23+& 00201 (cosIVals(i)*cosJVals(j)+cosIVals(i)*cosKVals(k)+c23)*3.0_dp+& 00202 (cosIVals(i)+cosJVals(j)+cosKVals(k))*9.0_dp+& 00203 27.0_dp) 00204 00205 spline_g%cc(ii)=spline_g%cc(ii)*coeff 00206 00207 END DO 00208 DEALLOCATE(cosIVals, cosJVals, cosKVals, stat=stat) 00209 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00210 00211 CALL timestop(handle) 00212 END SUBROUTINE pw_spline2_interpolate_values_g 00213 00214 ! ***************************************************************************** 00229 SUBROUTINE pw_spline3_interpolate_values_g(spline_g,error) 00230 TYPE(pw_type), POINTER :: spline_g 00231 TYPE(cp_error_type), INTENT(inout) :: error 00232 00233 CHARACTER(len=*), PARAMETER :: 00234 routineN = 'pw_spline3_interpolate_values_g', 00235 routineP = moduleN//':'//routineN 00236 00237 INTEGER :: handle, i, ii, j, k, stat 00238 INTEGER, DIMENSION(2, 3) :: gbo 00239 INTEGER, DIMENSION(3) :: n_tot 00240 LOGICAL :: failure 00241 REAL(KIND=dp) :: c23, coeff 00242 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cosIVals, cosJVals, cosKVals 00243 00244 CALL timeset(routineN,handle) 00245 00246 failure=.FALSE. 00247 n_tot(1:3) = spline_g%pw_grid%npts (1:3) 00248 gbo = spline_g%pw_grid%bounds 00249 00250 CPPrecondition(ASSOCIATED(spline_g),cp_failure_level,routineP,error,failure) 00251 CPPrecondition(spline_g%in_use==COMPLEXDATA1D,cp_failure_level,routineP,error,failure) 00252 CPPrecondition(spline_g%in_space==RECIPROCALSPACE,cp_failure_level,routineP,error,failure) 00253 CPPrecondition(.NOT.spline_g%pw_grid%spherical,cp_failure_level,routineP,error,failure) 00254 CPPrecondition(spline_g%pw_grid%grid_span==FULLSPACE,cp_failure_level,routineP,error,failure) 00255 00256 ALLOCATE(cosIVals(gbo(1,1):gbo(2,1)),& 00257 cosJVals(gbo(1,2):gbo(2,2)),& 00258 cosKVals(gbo(1,3):gbo(2,3)),stat=stat) 00259 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00260 00261 coeff=twopi/n_tot(1) 00262 !$omp parallel do private(i) 00263 DO i=gbo(1,1),gbo(2,1) 00264 cosIVals(i)=COS(coeff*REAL(i,dp)) 00265 END DO 00266 coeff=twopi/n_tot(2) 00267 !$omp parallel do private(j) 00268 DO j=gbo(1,2),gbo(2,2) 00269 cosJVals(j)=COS(coeff*REAL(j,dp)) 00270 END DO 00271 coeff=twopi/n_tot(3) 00272 !$omp parallel do private(k) 00273 DO k=gbo(1,3),gbo(2,3) 00274 cosKVals(k)=COS(coeff*REAL(k,dp)) 00275 END DO 00276 00277 !$omp parallel do private(i,j,k,ii,coeff,c23) 00278 DO ii=1,SIZE(spline_g%cc) 00279 i=spline_g%pw_grid%g_hat(1,ii) 00280 j=spline_g%pw_grid%g_hat(2,ii) 00281 k=spline_g%pw_grid%g_hat(3,ii) 00282 ! no opt 00283 !FM coeff=1.0/((cosVal(1)*cosVal(2)*cosVal(3))/27.0_dp+& 00284 !FM (cosVal(1)*cosVal(2)+cosVal(1)*cosVal(3)+& 00285 !FM cosVal(2)*cosVal(3))*2.0_dp/27.0_dp+& 00286 !FM (cosVal(1)+cosVal(2)+cosVal(3))*4.0_dp/27.0_dp+& 00287 !FM 8.0_dp/27.0_dp) 00288 ! opt 00289 c23=cosJVals(j)*cosKVals(k) 00290 coeff=27.0_dp/(cosIVals(i)*c23+& 00291 (cosIVals(i)*cosJVals(j)+cosIVals(i)*cosKVals(k)+c23)*2.0_dp+& 00292 (cosIVals(i)+cosJVals(j)+cosKVals(k))*4.0_dp+& 00293 8.0_dp) 00294 00295 spline_g%cc(ii)=spline_g%cc(ii)*coeff 00296 00297 END DO 00298 DEALLOCATE(cosIVals, cosJVals, cosKVals, stat=stat) 00299 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00300 00301 CALL timestop(handle) 00302 END SUBROUTINE pw_spline3_interpolate_values_g 00303 00304 ! ***************************************************************************** 00316 SUBROUTINE pw_spline2_evaluate_values_g(spline_g,error) 00317 TYPE(pw_type), POINTER :: spline_g 00318 TYPE(cp_error_type), INTENT(inout) :: error 00319 00320 CHARACTER(len=*), PARAMETER :: routineN = 'pw_spline2_evaluate_values_g', 00321 routineP = moduleN//':'//routineN 00322 00323 INTEGER :: handle, i, ii, j, k, stat 00324 INTEGER, DIMENSION(2, 3) :: bo, gbo 00325 INTEGER, DIMENSION(3) :: n, n_tot 00326 LOGICAL :: failure 00327 REAL(KIND=dp) :: c23, coeff, inv64 00328 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cosIVals, cosJVals, cosKVals 00329 00330 CALL timeset(routineN,handle) 00331 00332 failure=.FALSE. 00333 n(1:3) = spline_g%pw_grid%npts_local (1:3) 00334 n_tot(1:3) = spline_g%pw_grid%npts (1:3) 00335 bo = spline_g%pw_grid%bounds_local 00336 gbo = spline_g%pw_grid%bounds 00337 inv64=1.0_dp/64.0_dp 00338 00339 CPPrecondition(ASSOCIATED(spline_g),cp_failure_level,routineP,error,failure) 00340 CPPrecondition(spline_g%in_use==COMPLEXDATA1D,cp_failure_level,routineP,error,failure) 00341 CPPrecondition(spline_g%in_space==RECIPROCALSPACE,cp_failure_level,routineP,error,failure) 00342 CPPrecondition(.NOT.spline_g%pw_grid%spherical,cp_failure_level,routineP,error,failure) 00343 CPPrecondition(spline_g%pw_grid%grid_span==FULLSPACE,cp_failure_level,routineP,error,failure) 00344 00345 ALLOCATE(cosIVals(gbo(1,1):gbo(2,1)),cosJVals(gbo(1,2):gbo(2,2)), & 00346 cosKVals(gbo(1,3):gbo(2,3)), stat=stat) 00347 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00348 00349 coeff=twopi/n_tot(1) 00350 !$omp parallel do private(i) 00351 DO i=gbo(1,1),gbo(2,1) 00352 cosIVals(i)=COS(coeff*REAL(i,dp)) 00353 END DO 00354 coeff=twopi/n_tot(2) 00355 !$omp parallel do private(j) 00356 DO j=gbo(1,2),gbo(2,2) 00357 cosJVals(j)=COS(coeff*REAL(j,dp)) 00358 END DO 00359 !$omp parallel do private(k) 00360 DO k=gbo(1,3),gbo(2,3) 00361 cosKVals(k)=COS(coeff*REAL(k,dp)) 00362 END DO 00363 00364 !$omp parallel do private(ii,k,j,i,coeff,c23) 00365 DO ii=1,SIZE(spline_g%cc) 00366 i=spline_g%pw_grid%g_hat(1,ii) 00367 j=spline_g%pw_grid%g_hat(2,ii) 00368 k=spline_g%pw_grid%g_hat(3,ii) 00369 00370 c23=cosJVals(j)*cosKVals(k) 00371 coeff=(cosIVals(i)*c23+& 00372 (cosIVals(i)*cosJVals(j)+cosIVals(i)*cosKVals(k)+c23)*3.0_dp+& 00373 (cosIVals(i)+cosJVals(j)+cosKVals(k))*9.0_dp+& 00374 27.0_dp)*inv64 00375 00376 spline_g%cc(ii)=spline_g%cc(ii)*coeff 00377 END DO 00378 DEALLOCATE(cosIVals, cosJVals, cosKVals, stat=stat) 00379 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00380 00381 CALL timestop(handle) 00382 END SUBROUTINE pw_spline2_evaluate_values_g 00383 00384 ! ***************************************************************************** 00404 SUBROUTINE pw_nn_smear_g(spline_g,coeffs,error) 00405 TYPE(pw_type), POINTER :: spline_g 00406 REAL(KIND=dp), DIMENSION(4), INTENT(in) :: coeffs 00407 TYPE(cp_error_type), INTENT(inout) :: error 00408 00409 CHARACTER(len=*), PARAMETER :: routineN = 'pw_nn_smear_g', 00410 routineP = moduleN//':'//routineN 00411 00412 INTEGER :: handle, i, ii, j, k, stat 00413 INTEGER, DIMENSION(2, 3) :: bo, gbo 00414 INTEGER, DIMENSION(3) :: n, n_tot 00415 LOGICAL :: failure 00416 REAL(KIND=dp) :: c23, coeff 00417 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cosIVals, cosJVals, cosKVals 00418 REAL(KIND=dp), DIMENSION(4) :: r_coeffs 00419 00420 CALL timeset(routineN,handle) 00421 00422 failure=.FALSE. 00423 n(1:3) = spline_g%pw_grid%npts_local (1:3) 00424 n_tot(1:3) = spline_g%pw_grid%npts (1:3) 00425 bo = spline_g%pw_grid%bounds_local 00426 gbo = spline_g%pw_grid%bounds 00427 00428 r_coeffs=coeffs 00429 r_coeffs(2)=r_coeffs(2)*2.0_dp 00430 r_coeffs(3)=r_coeffs(3)*4.0_dp 00431 r_coeffs(4)=r_coeffs(4)*8.0_dp 00432 00433 CPPrecondition(ASSOCIATED(spline_g),cp_failure_level,routineP,error,failure) 00434 CPPrecondition(spline_g%in_use==COMPLEXDATA1D,cp_failure_level,routineP,error,failure) 00435 CPPrecondition(spline_g%in_space==RECIPROCALSPACE,cp_failure_level,routineP,error,failure) 00436 CPPrecondition(.NOT.spline_g%pw_grid%spherical,cp_failure_level,routineP,error,failure) 00437 CPPrecondition(spline_g%pw_grid%grid_span==FULLSPACE,cp_failure_level,routineP,error,failure) 00438 00439 ALLOCATE(cosIVals(gbo(1,1):gbo(2,1)),cosJVals(gbo(1,2):gbo(2,2)), & 00440 cosKVals(gbo(1,3):gbo(2,3)), stat=stat) 00441 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00442 00443 coeff=twopi/n_tot(1) 00444 !$omp parallel do private(i) 00445 DO i=gbo(1,1),gbo(2,1) 00446 cosIVals(i)=COS(coeff*REAL(i,dp)) 00447 END DO 00448 coeff=twopi/n_tot(2) 00449 !$omp parallel do private(j) 00450 DO j=gbo(1,2),gbo(2,2) 00451 cosJVals(j)=COS(coeff*REAL(j,dp)) 00452 END DO 00453 !$omp parallel do private(k) 00454 DO k=gbo(1,3),gbo(2,3) 00455 cosKVals(k)=COS(coeff*REAL(k,dp)) 00456 END DO 00457 00458 !$omp parallel do private(ii,k,j,i,coeff,c23) 00459 DO ii=1,SIZE(spline_g%cc) 00460 i=spline_g%pw_grid%g_hat(1,ii) 00461 j=spline_g%pw_grid%g_hat(2,ii) 00462 k=spline_g%pw_grid%g_hat(3,ii) 00463 00464 c23=cosJVals(j)*cosKVals(k) 00465 coeff=(r_coeffs(4)*cosIVals(i)*c23+& 00466 (cosIVals(i)*cosJVals(j)+cosIVals(i)*cosKVals(k)+c23)*r_coeffs(3)+& 00467 (cosIVals(i)+cosJVals(j)+cosKVals(k))*r_coeffs(2)+& 00468 r_coeffs(1)) 00469 00470 spline_g%cc(ii)=spline_g%cc(ii)*coeff 00471 END DO 00472 DEALLOCATE(cosIVals, cosJVals, cosKVals, stat=stat) 00473 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00474 00475 CALL timestop(handle) 00476 END SUBROUTINE pw_nn_smear_g 00477 00478 ! ***************************************************************************** 00490 SUBROUTINE pw_spline3_evaluate_values_g(spline_g,error) 00491 TYPE(pw_type), POINTER :: spline_g 00492 TYPE(cp_error_type), INTENT(inout) :: error 00493 00494 CHARACTER(len=*), PARAMETER :: routineN = 'pw_spline3_evaluate_values_g', 00495 routineP = moduleN//':'//routineN 00496 00497 INTEGER :: handle, i, ii, j, k, stat 00498 INTEGER, DIMENSION(2, 3) :: bo, gbo 00499 INTEGER, DIMENSION(3) :: n, n_tot 00500 LOGICAL :: failure 00501 REAL(KIND=dp) :: c23, coeff, inv27 00502 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cosIVals, cosJVals, cosKVals 00503 00504 CALL timeset(routineN,handle) 00505 00506 failure=.FALSE. 00507 n(1:3) = spline_g%pw_grid%npts_local (1:3) 00508 n_tot(1:3) = spline_g%pw_grid%npts (1:3) 00509 bo = spline_g%pw_grid%bounds_local 00510 gbo = spline_g%pw_grid%bounds 00511 inv27=1.0_dp/27.0_dp 00512 00513 CPPrecondition(ASSOCIATED(spline_g),cp_failure_level,routineP,error,failure) 00514 CPPrecondition(spline_g%in_use==COMPLEXDATA1D,cp_failure_level,routineP,error,failure) 00515 CPPrecondition(spline_g%in_space==RECIPROCALSPACE,cp_failure_level,routineP,error,failure) 00516 CPPrecondition(.NOT.spline_g%pw_grid%spherical,cp_failure_level,routineP,error,failure) 00517 CPPrecondition(spline_g%pw_grid%grid_span==FULLSPACE,cp_failure_level,routineP,error,failure) 00518 00519 ALLOCATE(cosIVals(gbo(1,1):gbo(2,1)),cosJVals(gbo(1,2):gbo(2,2)), & 00520 cosKVals(gbo(1,3):gbo(2,3)), stat=stat) 00521 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00522 00523 coeff=twopi/n_tot(1) 00524 !$omp parallel do private(i) 00525 DO i=gbo(1,1),gbo(2,1) 00526 cosIVals(i)=COS(coeff*REAL(i,dp)) 00527 END DO 00528 coeff=twopi/n_tot(2) 00529 !$omp parallel do private(j) 00530 DO j=gbo(1,2),gbo(2,2) 00531 cosJVals(j)=COS(coeff*REAL(j,dp)) 00532 END DO 00533 coeff=twopi/n_tot(3) 00534 !$omp parallel do private(k) 00535 DO k=gbo(1,3),gbo(2,3) 00536 cosKVals(k)=COS(coeff*REAL(k,dp)) 00537 END DO 00538 00539 !$omp parallel do private(ii,k,j,i,coeff,c23) 00540 DO ii=1,SIZE(spline_g%cc) 00541 i=spline_g%pw_grid%g_hat(1,ii) 00542 j=spline_g%pw_grid%g_hat(2,ii) 00543 k=spline_g%pw_grid%g_hat(3,ii) 00544 ! no opt 00545 !FM coeff=((cosVal(1)*cosVal(2)*cosVal(3))/27.0_dp+& 00546 !FM (cosVal(1)*cosVal(2)+cosVal(1)*cosVal(3)+& 00547 !FM cosVal(2)*cosVal(3))*2.0_dp/27.0_dp+& 00548 !FM (cosVal(1)+cosVal(2)+cosVal(3))*4.0_dp/27.0_dp+& 00549 !FM 8.0_dp/27.0_dp) 00550 ! opt 00551 c23=cosJVals(j)*cosKVals(k) 00552 coeff=(cosIVals(i)*c23+& 00553 (cosIVals(i)*cosJVals(j)+cosIVals(i)*cosKVals(k)+c23)*2.0_dp+& 00554 (cosIVals(i)+cosJVals(j)+cosKVals(k))*4.0_dp+& 00555 8.0_dp)*inv27 00556 00557 spline_g%cc(ii)=spline_g%cc(ii)*coeff 00558 END DO 00559 DEALLOCATE(cosIVals, cosJVals, cosKVals, stat=stat) 00560 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00561 00562 CALL timestop(handle) 00563 END SUBROUTINE pw_spline3_evaluate_values_g 00564 00565 ! ***************************************************************************** 00578 SUBROUTINE pw_spline_scale_deriv(deriv_vals_r,cell,transpose,scale,error) 00579 TYPE(pw_p_type), DIMENSION(3) :: deriv_vals_r 00580 TYPE(cell_type), POINTER :: cell 00581 LOGICAL, INTENT(in), OPTIONAL :: transpose 00582 REAL(KIND=dp), INTENT(in), OPTIONAL :: scale 00583 TYPE(cp_error_type), INTENT(inout) :: error 00584 00585 CHARACTER(len=*), PARAMETER :: routineN = 'pw_spline_scale_deriv', 00586 routineP = moduleN//':'//routineN 00587 00588 INTEGER :: handle, i, idir, j, k 00589 INTEGER, DIMENSION(2, 3) :: bo 00590 INTEGER, DIMENSION(3) :: n_tot 00591 LOGICAL :: diag, failure, my_transpose 00592 REAL(KIND=dp) :: dVal1, dVal2, dVal3, 00593 my_scale, scalef 00594 REAL(KIND=dp), DIMENSION(3, 3) :: h_grid 00595 REAL(kind=dp), DIMENSION(:, :, :), 00596 POINTER :: ddata, ddata2, ddata3 00597 00598 ! sun does not like that its components are modified 00599 00600 CALL timeset(routineN,handle) 00601 00602 failure=.FALSE. 00603 my_transpose=.FALSE. 00604 IF (PRESENT(transpose)) my_transpose=transpose 00605 my_scale=1.0_dp 00606 IF (PRESENT(scale)) my_scale=scale 00607 n_tot(1:3) = deriv_vals_r(1)%pw%pw_grid%npts (1:3) 00608 bo = deriv_vals_r(1)%pw%pw_grid%bounds_local 00609 00610 CPPrecondition(ASSOCIATED(cell),cp_failure_level,routineP,error,failure) 00611 00612 ! map grid to real derivative 00613 diag=.TRUE. 00614 IF (my_transpose) THEN 00615 DO j=1,3 00616 DO i=1,3 00617 h_grid(j,i)=my_scale*REAL(n_tot(i),dp)*cell%h_inv(i,j) 00618 IF (i/=j.AND.h_grid(j,i)/=0.0_dp) diag=.FALSE. 00619 END DO 00620 END DO 00621 ELSE 00622 DO j=1,3 00623 DO i=1,3 00624 h_grid(i,j)=my_scale*REAL(n_tot(i),dp)*cell%h_inv(i,j) 00625 IF (i/=j.AND.h_grid(i,j)/=0.0_dp) diag=.FALSE. 00626 END DO 00627 END DO 00628 END IF 00629 00630 IF (diag) THEN 00631 DO idir=1,3 00632 ddata => deriv_vals_r(idir)%pw%cr3d 00633 scalef=h_grid(idir,idir) 00634 CALL dscal((bo(2,1)-bo(1,1)+1)*(bo(2,2)-bo(1,2)+1)*(bo(2,3)-bo(1,3)+1),& 00635 scalef,ddata,1) 00636 !FM !$omp parallel do default(none) private(k,j,i) shared(ddata,scalef) 00637 !FM DO k = bo(1,3), bo(2,3) 00638 !FM DO j = bo(1,2), bo(2,2) 00639 !FM DO i = bo(1,1), bo(2,1) 00640 !FM ddata(i,j,k)=scalef*ddata(i,j,k) 00641 !FM END DO 00642 !FM END DO 00643 !FM END DO 00644 END DO 00645 ELSE 00646 ddata => deriv_vals_r(1)%pw%cr3d 00647 ddata2 => deriv_vals_r(2)%pw%cr3d 00648 ddata3 => deriv_vals_r(3)%pw%cr3d 00649 !$omp parallel do default(none) private(k,j,i,dVal1,dVal2,dVal3) & 00650 !$omp shared(ddata,ddata2,ddata3,h_grid,bo) 00651 DO k = bo(1,3), bo(2,3) 00652 DO j = bo(1,2), bo(2,2) 00653 DO i = bo(1,1), bo(2,1) 00654 00655 dVal1=ddata(i,j,k) 00656 dVal2=ddata2(i,j,k) 00657 dVal3=ddata3(i,j,k) 00658 00659 ddata(i,j,k)=h_grid(1,1)*dVal1+& 00660 h_grid(2,1)*dVal2+ h_grid(3,1)*dVal3 00661 ddata2(i,j,k)=h_grid(1,2)*dVal1+& 00662 h_grid(2,2)*dVal2+ h_grid(3,2)*dVal3 00663 ddata3(i,j,k)=h_grid(1,3)*dVal1+& 00664 h_grid(2,3)*dVal2+ h_grid(3,3)*dVal3 00665 00666 END DO 00667 END DO 00668 END DO 00669 END IF 00670 00671 CALL timestop(handle) 00672 END SUBROUTINE pw_spline_scale_deriv 00673 00674 ! ***************************************************************************** 00688 SUBROUTINE pw_spline3_deriv_g(spline_g,idir,error) 00689 TYPE(pw_type), POINTER :: spline_g 00690 INTEGER, INTENT(in) :: idir 00691 TYPE(cp_error_type), INTENT(inout) :: error 00692 00693 CHARACTER(len=*), PARAMETER :: routineN = 'pw_spline3_deriv_g', 00694 routineP = moduleN//':'//routineN 00695 00696 INTEGER :: handle, i, ii, j, k, stat 00697 INTEGER, DIMENSION(2, 3) :: bo, gbo 00698 INTEGER, DIMENSION(3) :: n, n_tot 00699 LOGICAL :: failure 00700 REAL(KIND=dp) :: coeff, inv9, tmp 00701 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: csIVals, csJVals, csKVals 00702 00703 CALL timeset(routineN,handle) 00704 00705 failure=.FALSE. 00706 n(1:3) = spline_g%pw_grid%npts_local (1:3) 00707 n_tot(1:3) = spline_g%pw_grid%npts (1:3) 00708 bo = spline_g%pw_grid%bounds_local 00709 gbo = spline_g%pw_grid%bounds 00710 inv9=1.0_dp/9.0_dp 00711 00712 CPPrecondition(ASSOCIATED(spline_g),cp_failure_level,routineP,error,failure) 00713 CPPrecondition(spline_g%in_use==COMPLEXDATA1D,cp_failure_level,routineP,error,failure) 00714 CPPrecondition(spline_g%in_space==RECIPROCALSPACE,cp_failure_level,routineP,error,failure) 00715 CPPrecondition(.not.spline_g%pw_grid%spherical,cp_failure_level,routineP,error,failure) 00716 CPPrecondition(spline_g%pw_grid%grid_span==FULLSPACE,cp_failure_level,routineP,error,failure) 00717 00718 ALLOCATE(csIVals(gbo(1,1):gbo(2,1)),& 00719 csJVals(gbo(1,2):gbo(2,2)),& 00720 csKVals(gbo(1,3):gbo(2,3)), stat=stat) 00721 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00722 00723 coeff=twopi/n_tot(1) 00724 IF (idir==1) THEN 00725 !$omp parallel do private(i) 00726 DO i=gbo(1,1),gbo(2,1) 00727 csIVals(i)=SIN(coeff*REAL(i,dp)) 00728 END DO 00729 ELSE 00730 !$omp parallel do private(i) 00731 DO i=gbo(1,1),gbo(2,1) 00732 csIVals(i)=COS(coeff*REAL(i,dp)) 00733 END DO 00734 END IF 00735 coeff=twopi/n_tot(2) 00736 IF (idir==2) THEN 00737 !$omp parallel do private(j) 00738 DO j=gbo(1,2),gbo(2,2) 00739 csJVals(j)=SIN(coeff*REAL(j,dp)) 00740 END DO 00741 ELSE 00742 !$omp parallel do private(j) 00743 DO j=gbo(1,2),gbo(2,2) 00744 csJVals(j)=COS(coeff*REAL(j,dp)) 00745 END DO 00746 END IF 00747 coeff=twopi/n_tot(3) 00748 IF (idir==3) THEN 00749 !$omp parallel do private(k) 00750 DO k=gbo(1,3),gbo(2,3) 00751 csKVals(k)=SIN(coeff*REAL(k,dp)) 00752 END DO 00753 ELSE 00754 !$omp parallel do private(k) 00755 DO k=gbo(1,3),gbo(2,3) 00756 csKVals(k)=COS(coeff*REAL(k,dp)) 00757 END DO 00758 END IF 00759 00760 SELECT CASE(idir) 00761 CASE (1) 00762 ! x deriv 00763 !$omp parallel do private(ii,k,j,i,coeff,tmp) 00764 DO ii=1,SIZE(spline_g%cc) 00765 i=spline_g%pw_grid%g_hat(1,ii) 00766 j=spline_g%pw_grid%g_hat(2,ii) 00767 k=spline_g%pw_grid%g_hat(3,ii) 00768 !FM ! formula 00769 !FM coeff=(sinVal(1)*cosVal(2)*cosVal(3))/9.0_dp+& 00770 !FM (sinVal(1)*cosVal(2)+sinVal(1)*cosVal(3))*2.0_dp/9.0_dp+& 00771 !FM sinVal(1)*4.0_dp/9.0_dp 00772 tmp=csIVals(i)*csJVals(j) 00773 coeff=(tmp*csKVals(k)+& 00774 (tmp+csIVals(i)*csKVals(k))*2.0_dp+& 00775 csIVals(i)*4.0_dp)*inv9 00776 00777 spline_g%cc(ii)=spline_g%cc(ii)*& 00778 CMPLX(0.0_dp,coeff,dp) 00779 END DO 00780 CASE(2) 00781 ! y deriv 00782 !$omp parallel do private(ii,k,j,i,coeff,tmp) 00783 DO ii=1,SIZE(spline_g%cc) 00784 i=spline_g%pw_grid%g_hat(1,ii) 00785 j=spline_g%pw_grid%g_hat(2,ii) 00786 k=spline_g%pw_grid%g_hat(3,ii) 00787 00788 tmp=csIVals(i)*csJVals(j) 00789 coeff=(tmp*csKVals(k)+& 00790 (tmp+csJVals(j)*csKVals(k))*2.0_dp+& 00791 csJVals(j)*4.0_dp)*inv9 00792 00793 spline_g%cc(ii)=spline_g%cc(ii)*& 00794 CMPLX(0.0_dp,coeff,dp) 00795 END DO 00796 CASE(3) 00797 ! z deriv 00798 !$omp parallel do private(ii,k,j,i,coeff,tmp) 00799 DO ii=1,SIZE(spline_g%cc) 00800 i=spline_g%pw_grid%g_hat(1,ii) 00801 j=spline_g%pw_grid%g_hat(2,ii) 00802 k=spline_g%pw_grid%g_hat(3,ii) 00803 00804 tmp=csIVals(i)*csKVals(k) 00805 coeff=(tmp*csJVals(j)+& 00806 (tmp+csJVals(j)*csKVals(k))*2.0_dp+& 00807 csKVals(k)*4.0_dp)*inv9 00808 00809 spline_g%cc(ii)=spline_g%cc(ii)*& 00810 CMPLX(0.0_dp,coeff,dp) 00811 END DO 00812 END SELECT 00813 00814 DEALLOCATE(csIVals, csJVals, csKVals, stat=stat) 00815 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00816 00817 CALL timestop(handle) 00818 END SUBROUTINE pw_spline3_deriv_g 00819 00820 ! ***************************************************************************** 00834 SUBROUTINE pw_spline2_deriv_g(spline_g,idir,error) 00835 TYPE(pw_type), POINTER :: spline_g 00836 INTEGER, INTENT(in) :: idir 00837 TYPE(cp_error_type), INTENT(inout) :: error 00838 00839 CHARACTER(len=*), PARAMETER :: routineN = 'pw_spline2_deriv_g', 00840 routineP = moduleN//':'//routineN 00841 00842 INTEGER :: handle, i, ii, j, k, stat 00843 INTEGER, DIMENSION(2, 3) :: bo 00844 INTEGER, DIMENSION(3) :: n, n_tot 00845 LOGICAL :: failure 00846 REAL(KIND=dp) :: coeff, inv16, tmp 00847 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: csIVals, csJVals, csKVals 00848 00849 CALL timeset(routineN,handle) 00850 00851 failure=.FALSE. 00852 n(1:3) = spline_g%pw_grid%npts_local (1:3) 00853 n_tot(1:3) = spline_g%pw_grid%npts (1:3) 00854 bo = spline_g%pw_grid%bounds 00855 inv16=1.0_dp/16.0_dp 00856 00857 CPPrecondition(ASSOCIATED(spline_g),cp_failure_level,routineP,error,failure) 00858 CPPrecondition(spline_g%in_use==COMPLEXDATA1D,cp_failure_level,routineP,error,failure) 00859 CPPrecondition(spline_g%in_space==RECIPROCALSPACE,cp_failure_level,routineP,error,failure) 00860 CPPrecondition(.not.spline_g%pw_grid%spherical,cp_failure_level,routineP,error,failure) 00861 CPPrecondition(spline_g%pw_grid%grid_span==FULLSPACE,cp_failure_level,routineP,error,failure) 00862 00863 ALLOCATE(csIVals(bo(1,1):bo(2,1)),csJVals(bo(1,2):bo(2,2)),& 00864 csKVals(bo(1,3):bo(2,3)), stat=stat) 00865 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00866 00867 coeff=twopi/n_tot(1) 00868 IF (idir==1) THEN 00869 !$omp parallel do private(i) 00870 DO i=bo(1,1),bo(2,1) 00871 csIVals(i)=SIN(coeff*REAL(i,dp)) 00872 END DO 00873 ELSE 00874 !$omp parallel do private(i) 00875 DO i=bo(1,1),bo(2,1) 00876 csIVals(i)=COS(coeff*REAL(i,dp)) 00877 END DO 00878 END IF 00879 coeff=twopi/n_tot(2) 00880 IF (idir==2) THEN 00881 !$omp parallel do private(j) 00882 DO j=bo(1,2),bo(2,2) 00883 csJVals(j)=SIN(coeff*REAL(j,dp)) 00884 END DO 00885 ELSE 00886 !$omp parallel do private(j) 00887 DO j=bo(1,2),bo(2,2) 00888 csJVals(j)=COS(coeff*REAL(j,dp)) 00889 END DO 00890 END IF 00891 coeff=twopi/n_tot(3) 00892 IF (idir==3) THEN 00893 !$omp parallel do private(k) 00894 DO k=bo(1,3),bo(2,3) 00895 csKVals(k)=SIN(coeff*REAL(k,dp)) 00896 END DO 00897 ELSE 00898 !$omp parallel do private(k) 00899 DO k=bo(1,3),bo(2,3) 00900 csKVals(k)=COS(coeff*REAL(k,dp)) 00901 END DO 00902 END IF 00903 00904 SELECT CASE(idir) 00905 CASE (1) 00906 ! x deriv 00907 !$omp parallel do private(ii,k,j,i,coeff,tmp) 00908 DO ii=1,SIZE(spline_g%cc) 00909 i=spline_g%pw_grid%g_hat(1,ii) 00910 j=spline_g%pw_grid%g_hat(2,ii) 00911 k=spline_g%pw_grid%g_hat(3,ii) 00912 !FM ! formula 00913 !FM coeff=(sinVal(1)*cosVal(2)*cosVal(3))/16.0_dp+& 00914 !FM (sinVal(1)*cosVal(2)+sinVal(1)*cosVal(3))*3.0_dp/16.0_dp+& 00915 !FM sinVal(1)*9.0_dp/16.0_dp 00916 tmp=csIVals(i)*csJVals(j) 00917 coeff=(tmp*csKVals(k)+& 00918 (tmp+csIVals(i)*csKVals(k))*3.0_dp+& 00919 csIVals(i)*9.0_dp)*inv16 00920 00921 spline_g%cc(ii)=spline_g%cc(ii)*& 00922 CMPLX(0.0_dp,coeff,dp) 00923 END DO 00924 CASE(2) 00925 ! y deriv 00926 !$omp parallel do private(ii,k,j,i,coeff,tmp) 00927 DO ii=1,SIZE(spline_g%cc) 00928 i=spline_g%pw_grid%g_hat(1,ii) 00929 j=spline_g%pw_grid%g_hat(2,ii) 00930 k=spline_g%pw_grid%g_hat(3,ii) 00931 00932 tmp=csIVals(i)*csJVals(j) 00933 coeff=(tmp*csKVals(k)+& 00934 (tmp+csJVals(j)*csKVals(k))*3.0_dp+& 00935 csJVals(j)*9.0_dp)*inv16 00936 00937 spline_g%cc(ii)=spline_g%cc(ii)*& 00938 CMPLX(0.0_dp,coeff,dp) 00939 END DO 00940 CASE(3) 00941 ! z deriv 00942 !$omp parallel do private(ii,k,j,i,coeff,tmp) 00943 DO ii=1,SIZE(spline_g%cc) 00944 i=spline_g%pw_grid%g_hat(1,ii) 00945 j=spline_g%pw_grid%g_hat(2,ii) 00946 k=spline_g%pw_grid%g_hat(3,ii) 00947 00948 tmp=csIVals(i)*csKVals(k) 00949 coeff=(tmp*csJVals(j)+& 00950 (tmp+csJVals(j)*csKVals(k))*3.0_dp+& 00951 csKVals(k)*9.0_dp)*inv16 00952 00953 spline_g%cc(ii)=spline_g%cc(ii)*& 00954 CMPLX(0.0_dp,coeff,dp) 00955 END DO 00956 END SELECT 00957 00958 DEALLOCATE(csIVals, csJVals, csKVals, stat=stat) 00959 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00960 00961 CALL timestop(handle) 00962 END SUBROUTINE pw_spline2_deriv_g 00963 00964 ! ***************************************************************************** 00981 SUBROUTINE pw_compose_stripe(weights,in_val,in_val_first,in_val_last,& 00982 out_val,n_el) 00983 REAL(kind=dp), DIMENSION(0:2), 00984 INTENT(in) :: weights 00985 REAL(kind=dp), DIMENSION(*), INTENT(in) :: in_val 00986 REAL(kind=dp), INTENT(in) :: in_val_first, in_val_last 00987 REAL(kind=dp), DIMENSION(*), 00988 INTENT(inout) :: out_val 00989 INTEGER :: n_el 00990 00991 CHARACTER(len=*), PARAMETER :: routineN = 'pw_compose_stripe', 00992 routineP = moduleN//':'//routineN 00993 00994 INTEGER :: i 00995 LOGICAL :: failure 00996 REAL(kind=dp) :: v0, v1, v2 00997 00998 !1:n_el), & 00999 !1:n_el), & 01000 01001 IF (n_el<1) RETURN 01002 failure=.FALSE. 01003 v0=in_val_first 01004 v1=in_val(1) 01005 IF (weights(1)==0.0_dp) THEN 01006 ! optimized version for x deriv 01007 DO i=1,n_el-3,3 01008 v2=in_val(i+1) 01009 out_val(i)= out_val(i)+& 01010 weights(0)*v0+& 01011 weights(2)*v2 01012 v0=in_val(i+2) 01013 out_val(i+1)= out_val(i+1)+& 01014 weights(0)*v1+& 01015 weights(2)*v0 01016 v1=in_val(i+3) 01017 out_val(i+2)= out_val(i+2)+& 01018 weights(0)*v2+& 01019 weights(2)*v1 01020 END DO 01021 ELSE 01022 ! generic version 01023 DO i=1,n_el-3,3 01024 v2=in_val(i+1) 01025 out_val(i)= out_val(i)+& 01026 weights(0)*v0+& 01027 weights(1)*v1+& 01028 weights(2)*v2 01029 v0=in_val(i+2) 01030 out_val(i+1)= out_val(i+1)+& 01031 weights(0)*v1+& 01032 weights(1)*v2+& 01033 weights(2)*v0 01034 v1=in_val(i+3) 01035 out_val(i+2)= out_val(i+2)+& 01036 weights(0)*v2+& 01037 weights(1)*v0+& 01038 weights(2)*v1 01039 END DO 01040 END IF 01041 SELECT CASE(MODULO(n_el-1,3)) 01042 CASE(0) 01043 v2=in_val_last 01044 out_val(n_el)= out_val(n_el)+& 01045 weights(0)*v0+& 01046 weights(1)*v1+& 01047 weights(2)*v2 01048 CASE(1) 01049 v2=in_val(n_el) 01050 out_val(n_el-1)= out_val(n_el-1)+& 01051 weights(0)*v0+& 01052 weights(1)*v1+& 01053 weights(2)*v2 01054 v0=in_val_last 01055 out_val(n_el)= out_val(n_el)+& 01056 weights(0)*v1+& 01057 weights(1)*v2+& 01058 weights(2)*v0 01059 CASE(2) 01060 v2=in_val(n_el-1) 01061 out_val(n_el-2)= out_val(n_el-2)+& 01062 weights(0)*v0+& 01063 weights(1)*v1+& 01064 weights(2)*v2 01065 v0=in_val(n_el) 01066 out_val(n_el-1)= out_val(n_el-1)+& 01067 weights(0)*v1+& 01068 weights(1)*v2+& 01069 weights(2)*v0 01070 v1=in_val_last 01071 out_val(n_el)= out_val(n_el)+& 01072 weights(0)*v2+& 01073 weights(1)*v0+& 01074 weights(2)*v1 01075 END SELECT 01076 01077 END SUBROUTINE pw_compose_stripe 01078 ! ***************************************************************************** 01079 SUBROUTINE pw_compose_stripe2(weights,in_val,in_val_first,in_val_last,& 01080 out_val,first_val,last_val,myj,myk,j,k) 01081 REAL(kind=dp), DIMENSION(0:2), 01082 INTENT(in) :: weights 01083 REAL(kind=dp), DIMENSION(:, :, :), 01084 POINTER :: in_val 01085 REAL(kind=dp), INTENT(in) :: in_val_first, in_val_last 01086 REAL(kind=dp), DIMENSION(:, :, :), 01087 POINTER :: out_val 01088 INTEGER, INTENT(in) :: first_val, last_val, myj, 01089 myk, j, k 01090 01091 CHARACTER(len=*), PARAMETER :: routineN = 'pw_compose_stripe2', 01092 routineP = moduleN//':'//routineN 01093 01094 INTEGER :: i, n_el 01095 LOGICAL :: failure 01096 REAL(kind=dp) :: v0, v1, v2 01097 01098 IF (last_val-first_val<0) RETURN 01099 failure=.FALSE. 01100 v0=in_val_first 01101 v1=in_val(first_val,myj,myk) 01102 IF (weights(1)==0.0_dp) THEN 01103 ! optimized version for x deriv 01104 DO i=first_val,last_val-3,3 01105 v2=in_val(i+1,myj,myk) 01106 out_val(i,j,k)= out_val(i,j,k)+& 01107 weights(0)*v0+& 01108 weights(2)*v2 01109 v0=in_val(i+2,myj,myk) 01110 out_val(i+1,j,k)= out_val(i+1,j,k)+& 01111 weights(0)*v1+& 01112 weights(2)*v0 01113 v1=in_val(i+3,myj,myk) 01114 out_val(i+2,j,k)= out_val(i+2,j,k)+& 01115 weights(0)*v2+& 01116 weights(2)*v1 01117 END DO 01118 ELSE 01119 ! generic version 01120 DO i=first_val,last_val-3,3 01121 v2=in_val(i+1,myj,myk) 01122 out_val(i,j,k)= out_val(i,j,k)+& 01123 weights(0)*v0+& 01124 weights(1)*v1+& 01125 weights(2)*v2 01126 v0=in_val(i+2,myj,myk) 01127 out_val(i+1,j,k)= out_val(i+1,j,k)+& 01128 weights(0)*v1+& 01129 weights(1)*v2+& 01130 weights(2)*v0 01131 v1=in_val(i+3,myj,myk) 01132 out_val(i+2,j,k)= out_val(i+2,j,k)+& 01133 weights(0)*v2+& 01134 weights(1)*v0+& 01135 weights(2)*v1 01136 END DO 01137 END IF 01138 n_el=last_val 01139 SELECT CASE(MODULO(last_val-first_val,3)) 01140 CASE(0) 01141 v2=in_val_last 01142 out_val(n_el,j,k)= out_val(n_el,j,k)+& 01143 weights(0)*v0+& 01144 weights(1)*v1+& 01145 weights(2)*v2 01146 CASE(1) 01147 v2=in_val(n_el,myj,myk) 01148 out_val(n_el-1,j,k)= out_val(n_el-1,j,k)+& 01149 weights(0)*v0+& 01150 weights(1)*v1+& 01151 weights(2)*v2 01152 v0=in_val_last 01153 out_val(n_el,j,k)= out_val(n_el,j,k)+& 01154 weights(0)*v1+& 01155 weights(1)*v2+& 01156 weights(2)*v0 01157 CASE(2) 01158 v2=in_val(n_el-1,myj,myk) 01159 out_val(n_el-2,j,k)= out_val(n_el-2,j,k)+& 01160 weights(0)*v0+& 01161 weights(1)*v1+& 01162 weights(2)*v2 01163 v0=in_val(n_el,myj,myk) 01164 out_val(n_el-1,j,k)= out_val(n_el-1,j,k)+& 01165 weights(0)*v1+& 01166 weights(1)*v2+& 01167 weights(2)*v0 01168 v1=in_val_last 01169 out_val(n_el,j,k)= out_val(n_el,j,k)+& 01170 weights(0)*v2+& 01171 weights(1)*v0+& 01172 weights(2)*v1 01173 END SELECT 01174 01175 END SUBROUTINE pw_compose_stripe2 01176 !********** 01177 01178 ! ***************************************************************************** 01192 SUBROUTINE pw_nn_compose_r_work(weights,in_val,out_val,pw_in,bo,error) 01193 REAL(kind=dp), DIMENSION(0:2, 0:2, 0:2) :: weights 01194 INTEGER, DIMENSION(2, 3) :: bo 01195 TYPE(pw_type), POINTER :: pw_in 01196 REAL(kind=dp), DIMENSION(bo(1, 1):bo(2, & 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, & 3)), INTENT(inout) :: out_val 01197 REAL(kind=dp), DIMENSION(bo(1, 1):bo(2, & 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, & 3)), INTENT(in) :: in_val 01198 TYPE(cp_error_type), INTENT(inout) :: error 01199 01200 CHARACTER(len=*), PARAMETER :: routineN = 'pw_nn_compose_r_work', 01201 routineP = moduleN//':'//routineN 01202 01203 INTEGER :: i, j, jw, k, kw, myj, myk, 01204 stat 01205 INTEGER, DIMENSION(2, 3) :: gbo 01206 INTEGER, DIMENSION(3) :: s 01207 LOGICAL :: failure, has_boundary, 01208 yderiv, zderiv 01209 REAL(kind=dp) :: in_val_f, in_val_l 01210 REAL(kind=dp), DIMENSION(:, :), POINTER :: l_boundary, tmp, u_boundary 01211 01212 zderiv=ALL(weights(:,:,1)==0.0_dp) 01213 yderiv=ALL(weights(:,1,:)==0.0_dp) 01214 bo=pw_in%pw_grid%bounds_local 01215 gbo=pw_in%pw_grid%bounds 01216 DO i=1,3 01217 s(i)=bo(2,i)-bo(1,i)+1 01218 END DO 01219 IF (ANY(s<1)) RETURN 01220 has_boundary= ANY(pw_in%pw_grid%bounds_local(:,1) /= & 01221 pw_in%pw_grid%bounds(:,1)) 01222 IF (has_boundary) THEN 01223 ALLOCATE(l_boundary(bo(1,2):bo(2,2),bo(1,3):bo(2,3)),& 01224 u_boundary(bo(1,2):bo(2,2),bo(1,3):bo(2,3)),& 01225 tmp(bo(1,2):bo(2,2),bo(1,3):bo(2,3)),stat=stat) 01226 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01227 tmp(:,:)=pw_in%cr3d(bo(2,1),:,:) 01228 CALL mp_sendrecv(tmp,pw_in%pw_grid%para%pos_of_x(& 01229 gbo(1,1)+MODULO(bo(2,1)+1-gbo(1,1),gbo(2,1)-gbo(1,1)+1)),& 01230 l_boundary,pw_in%pw_grid%para%pos_of_x(& 01231 gbo(1,1)+MODULO(bo(1,1)-1-gbo(1,1),gbo(2,1)-gbo(1,1)+1)),& 01232 pw_in%pw_grid%para%group) 01233 tmp(:,:)=pw_in%cr3d(bo(1,1),:,:) 01234 CALL mp_sendrecv(tmp,pw_in%pw_grid%para%pos_of_x(& 01235 gbo(1,1)+MODULO(bo(1,1)-1-gbo(1,1),gbo(2,1)-gbo(1,1)+1)),& 01236 u_boundary,pw_in%pw_grid%para%pos_of_x(& 01237 gbo(1,1)+MODULO(bo(2,1)+1-gbo(1,1),gbo(2,1)-gbo(1,1)+1)),& 01238 pw_in%pw_grid%para%group) 01239 DEALLOCATE(tmp,stat=stat) 01240 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 01241 END IF 01242 01243 !$omp parallel do default(none) private(k,kw,myk,j,jw,myj,in_val_f,& 01244 !$omp in_val_l) shared(zderiv,yderiv,bo,in_val,out_val,s,l_boundary,& 01245 !$omp u_boundary,weights,has_boundary) 01246 DO k=0,s(3)-1 01247 DO kw=0,2 01248 myk=bo(1,3)+MODULO(k+kw-1,s(3)) 01249 IF (zderiv.AND.kw==1) CYCLE 01250 DO j=0,s(2)-1 01251 DO jw=0,2 01252 myj=bo(1,2)+MODULO(j+jw-1,s(2)) 01253 IF (yderiv.AND.jw==1) CYCLE 01254 IF (has_boundary) THEN 01255 in_val_f=l_boundary(myj,myk) 01256 in_val_l=u_boundary(myj,myk) 01257 ELSE 01258 in_val_f=in_val(bo(2,1),myj,myk) 01259 in_val_l=in_val(bo(1,1),myj,myk) 01260 END IF 01261 CALL pw_compose_stripe(weights=weights(:,jw,kw),& 01262 in_val=in_val(:,myj,myk),& 01263 in_val_first=in_val_f,in_val_last=in_val_l,& 01264 out_val=out_val(:,bo(1,2)+j,bo(1,3)+k),n_el=s(1)) 01265 END DO 01266 END DO 01267 END DO 01268 END DO 01269 IF (has_boundary) THEN 01270 DEALLOCATE(l_boundary,u_boundary,stat=stat) 01271 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 01272 END IF 01273 END SUBROUTINE pw_nn_compose_r_work 01274 01275 ! ***************************************************************************** 01286 SUBROUTINE pw_nn_compose_r(weights, pw_in, pw_out,error) 01287 REAL(kind=dp), DIMENSION(0:2, 0:2, 0:2) :: weights 01288 TYPE(pw_type), POINTER :: pw_in, pw_out 01289 TYPE(cp_error_type), INTENT(inout) :: error 01290 01291 CHARACTER(len=*), PARAMETER :: routineN = 'pw_nn_compose_r', 01292 routineP = moduleN//':'//routineN 01293 01294 INTEGER :: handle 01295 LOGICAL :: failure 01296 01297 CALL timeset(routineN,handle) 01298 failure=.FALSE. 01299 CPPrecondition(ASSOCIATED(pw_in),cp_failure_level,routineP,error,failure) 01300 CPPrecondition(pw_in%in_space==REALSPACE,cp_failure_level,routineP,error,failure) 01301 CPPrecondition(pw_in%in_use==REALDATA3D,cp_failure_level,routineP,error,failure) 01302 CPPrecondition(ASSOCIATED(pw_out),cp_failure_level,routineP,error,failure) 01303 CPPrecondition(pw_out%in_space==REALSPACE,cp_failure_level,routineP,error,failure) 01304 CPPrecondition(pw_out%in_use==REALDATA3D,cp_failure_level,routineP,error,failure) 01305 IF (.NOT.ALL(pw_in%pw_grid%bounds_local(:,2:3) == pw_in%pw_grid%bounds(:,2:3))) THEN 01306 CALL cp_assert(.FALSE.,& 01307 cp_failure_level,cp_assertion_failed,routineP,& 01308 "wrong pw distribution",error,failure) 01309 END IF 01310 IF (.NOT.failure) THEN 01311 CALL pw_nn_compose_r_work(weights=weights,in_val=pw_in%cr3d,& 01312 out_val=pw_out%cr3d,pw_in=pw_in,bo=pw_in%pw_grid%bounds_local,& 01313 error=error) 01314 END IF 01315 CALL timestop(handle) 01316 END SUBROUTINE pw_nn_compose_r 01317 01318 ! ***************************************************************************** 01336 SUBROUTINE pw_nn_smear_r(pw_in,pw_out,coeffs,error) 01337 TYPE(pw_type), POINTER :: pw_in, pw_out 01338 REAL(KIND=dp), DIMENSION(4), INTENT(in) :: coeffs 01339 TYPE(cp_error_type), INTENT(inout) :: error 01340 01341 CHARACTER(len=*), PARAMETER :: routineN = 'pw_nn_smear_r', 01342 routineP = moduleN//':'//routineN 01343 01344 INTEGER :: i, j, k 01345 REAL(kind=dp), 01346 DIMENSION(-1:1, -1:1, -1:1) :: weights 01347 01348 DO k=-1,1 01349 DO j=-1,1 01350 DO i=-1,1 01351 weights(i,j,k)=coeffs(ABS(i)+ABS(j)+ABS(k)+1) 01352 END DO 01353 END DO 01354 END DO 01355 01356 CALL pw_nn_compose_r(weights=weights,pw_in=pw_in,pw_out=pw_out,& 01357 error=error) 01358 END SUBROUTINE pw_nn_smear_r 01359 01360 ! ***************************************************************************** 01388 SUBROUTINE pw_nn_deriv_r(pw_in,pw_out,coeffs,idir,error) 01389 TYPE(pw_type), POINTER :: pw_in, pw_out 01390 REAL(KIND=dp), DIMENSION(3), INTENT(in) :: coeffs 01391 INTEGER :: idir 01392 TYPE(cp_error_type), INTENT(inout) :: error 01393 01394 CHARACTER(len=*), PARAMETER :: routineN = 'pw_nn_deriv_r', 01395 routineP = moduleN//':'//routineN 01396 01397 INTEGER :: i, idirVal, j, k 01398 LOGICAL :: failure 01399 REAL(kind=dp), 01400 DIMENSION(-1:1, -1:1, -1:1) :: weights 01401 01402 failure=.FALSE. 01403 DO k=-1,1 01404 DO j=-1,1 01405 DO i=-1,1 01406 SELECT CASE(idir) 01407 CASE (1) 01408 idirVal=i 01409 CASE (2) 01410 idirVal=j 01411 CASE (3) 01412 idirVal=k 01413 CASE default 01414 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 01415 routineP,"invalid idir ("//TRIM(cp_to_string(idir))//")",& 01416 error,failure) 01417 END SELECT 01418 IF (idirVal==0) THEN 01419 weights(i,j,k)=0.0_dp 01420 ELSE 01421 weights(i,j,k)=REAL(idirVal,dp)*coeffs(ABS(i)+ABS(j)+ABS(k)) 01422 END IF 01423 END DO 01424 END DO 01425 END DO 01426 01427 IF (.NOT.failure) THEN 01428 CALL pw_nn_compose_r(weights=weights,pw_in=pw_in,pw_out=pw_out,& 01429 error=error) 01430 END IF 01431 END SUBROUTINE pw_nn_deriv_r 01432 01433 ! ***************************************************************************** 01478 SUBROUTINE add_coarse2fine(coarse_coeffs_pw,fine_values_pw,& 01479 weights_1d,w_border0,w_border1,pbc,safe_computation,error) 01480 TYPE(pw_type), POINTER :: coarse_coeffs_pw, 01481 fine_values_pw 01482 REAL(kind=dp), DIMENSION(4), INTENT(in) :: weights_1d 01483 REAL(kind=dp), INTENT(in) :: w_border0 01484 REAL(kind=dp), DIMENSION(3), INTENT(in) :: w_border1 01485 LOGICAL, INTENT(in) :: pbc 01486 LOGICAL, INTENT(in), OPTIONAL :: safe_computation 01487 TYPE(cp_error_type), INTENT(inout) :: error 01488 01489 CHARACTER(len=*), PARAMETER :: routineN = 'add_coarse2fine', 01490 routineP = moduleN//':'//routineN 01491 01492 INTEGER :: coarse_slice_size, f_shift(3), fi, fi_lb, fi_ub, fj, fk, 01493 handle, handle2, i, ii, ij, ik, ip, j, k, my_lb, my_ub, n_procs, p, 01494 p_lb, p_old, p_ub, rcv_tot_size, rest_b, s(3), send_tot_size, sf, 01495 shift, ss, stat, x, x_att, xx 01496 INTEGER, ALLOCATABLE, DIMENSION(:) :: rcv_offset, rcv_size, 01497 real_rcv_size, send_offset, 01498 send_size, sent_size 01499 INTEGER, DIMENSION(2, 3) :: coarse_bo, coarse_gbo, 01500 fine_bo, fine_gbo, 01501 my_coarse_bo 01502 INTEGER, DIMENSION(:), POINTER :: pos_of_x 01503 LOGICAL :: failure, has_i_lbound, 01504 has_i_ubound, is_split, 01505 safe_calc 01506 REAL(kind=dp) :: v0, v1, v2, v3, wi, wj, wk 01507 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: rcv_buf, send_buf 01508 REAL(kind=dp), DIMENSION(3) :: w_0, ww0 01509 REAL(kind=dp), DIMENSION(4) :: w_1, ww1 01510 REAL(kind=dp), DIMENSION(:, :, :), 01511 POINTER :: coarse_coeffs, fine_values 01512 01513 CALL timeset(routineN,handle) 01514 ! CALL timeset(routineN//"_pre",handle2) 01515 failure=.FALSE. 01516 safe_calc=.FALSE. 01517 IF (PRESENT(safe_computation)) safe_calc=safe_computation 01518 CALL mp_comm_compare(coarse_coeffs_pw%pw_grid%para%group,& 01519 fine_values_pw%pw_grid%para%group,ii) 01520 IF (ii>1) THEN 01521 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 01522 END IF 01523 IF (.NOT. failure) THEN 01524 my_coarse_bo=coarse_coeffs_pw%pw_grid%bounds_local 01525 coarse_gbo=coarse_coeffs_pw%pw_grid%bounds 01526 fine_bo=fine_values_pw%pw_grid%bounds_local 01527 fine_gbo=fine_values_pw%pw_grid%bounds 01528 f_shift=fine_gbo(1,:)-2*coarse_gbo(1,:) 01529 DO j=2,3 01530 DO i=1,2 01531 coarse_bo(i,j)=FLOOR((fine_bo(i,j)-f_shift(j))/2.) 01532 END DO 01533 END DO 01534 IF (fine_bo(1,1)<=fine_bo(2,1)) THEN 01535 coarse_bo(1,1)=FLOOR((fine_bo(1,1)-2-f_shift(1))/2.) 01536 coarse_bo(2,1)=FLOOR((fine_bo(2,1)+3-f_shift(1))/2.) 01537 ELSE 01538 coarse_bo(1,1)=coarse_gbo(2,1) 01539 coarse_bo(2,1)=coarse_gbo(2,1)-1 01540 END IF 01541 is_split= ANY(coarse_gbo(:,1) /= my_coarse_bo(:,1)) 01542 IF (.NOT.is_split.OR..NOT.pbc) THEN 01543 coarse_bo(1,1)=MAX(coarse_gbo(1,1),coarse_bo(1,1)) 01544 coarse_bo(2,1)=MIN(coarse_gbo(2,1),coarse_bo(2,1)) 01545 END IF 01546 has_i_ubound=(fine_gbo(2,1)/=fine_bo(2,1)).or.pbc.and.is_split 01547 has_i_lbound=(fine_gbo(1,1)/=fine_bo(1,1)).or.pbc.and.is_split 01548 01549 IF (pbc) THEN 01550 CPPrecondition(ALL(fine_gbo(1,:)==2*coarse_gbo(1,:)+f_shift),cp_failure_level,routineP,error,failure) 01551 CPPrecondition(ALL(fine_gbo(2,:)==2*coarse_gbo(2,:)+1+f_shift),cp_failure_level,routineP,error,failure) 01552 ELSE 01553 CPPrecondition(ALL(fine_gbo(2,:)==2*coarse_gbo(2,:)+f_shift),cp_failure_level,routineP,error,failure) 01554 CPPrecondition(ALL(fine_gbo(1,:)==2*coarse_gbo(1,:)+f_shift),cp_failure_level,routineP,error,failure) 01555 END IF 01556 01557 coarse_coeffs => coarse_coeffs_pw%cr3d 01558 DO i=1,3 01559 s(i)=coarse_gbo(2,i)-coarse_gbo(1,i)+1 01560 END DO 01561 ! CALL timestop(handle2) 01562 ! *** parallel case 01563 IF (is_split) THEN 01564 CALL timeset(routineN//"_comm",handle2) 01565 coarse_slice_size=(coarse_bo(2,2)-coarse_bo(1,2)+1)*& 01566 (coarse_bo(2,3)-coarse_bo(1,3)+1) 01567 n_procs=coarse_coeffs_pw%pw_grid%para%group_size 01568 ALLOCATE(send_size(0:n_procs-1),send_offset(0:n_procs-1),& 01569 sent_size(0:n_procs-1),rcv_size(0:n_procs-1), & 01570 rcv_offset(0:n_procs-1),real_rcv_size(0:n_procs-1), stat=stat) 01571 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 01572 01573 ! ** rcv size count 01574 01575 pos_of_x => coarse_coeffs_pw%pw_grid%para%pos_of_x 01576 p_old=pos_of_x(coarse_gbo(1,1) & 01577 +MODULO(coarse_bo(1,1)-coarse_gbo(1,1),s(1))) 01578 rcv_size=0 01579 DO x=coarse_bo(1,1),coarse_bo(2,1) 01580 p=pos_of_x(coarse_gbo(1,1)+MODULO(x-coarse_gbo(1,1),s(1))) 01581 rcv_size(p)=rcv_size(p)+coarse_slice_size 01582 END DO 01583 01584 ! ** send size count 01585 01586 pos_of_x => fine_values_pw%pw_grid%para%pos_of_x 01587 sf=fine_gbo(2,1)-fine_gbo(1,1)+1 01588 fi_lb=2*my_coarse_bo(1,1)-3+f_shift(1) 01589 fi_ub=2*my_coarse_bo(2,1)+3+f_shift(1) 01590 IF (.not.pbc) THEN 01591 fi_lb=MAX(fi_lb,fine_gbo(1,1)) 01592 fi_ub=MIN(fi_ub,fine_gbo(2,1)) 01593 ELSE 01594 fi_ub=MIN(fi_ub,fi_lb+sf-1) 01595 END IF 01596 p_old=pos_of_x(fine_gbo(1,1)+MODULO(fi_lb-fine_gbo(1,1),sf)) 01597 p_lb=FLOOR((fi_lb-2-f_shift(1))/2.) 01598 send_size=0 01599 DO x=fi_lb,fi_ub 01600 p=pos_of_x(fine_gbo(1,1)+MODULO(x-fine_gbo(1,1),sf)) 01601 IF (p/=p_old) THEN 01602 p_ub=FLOOR((x-1+3-f_shift(1))/2.) 01603 01604 send_size(p_old)=send_size(p_old)+(MIN(p_ub,my_coarse_bo(2,1)) & 01605 -MAX(p_lb,my_coarse_bo(1,1))+1)*coarse_slice_size 01606 01607 IF (pbc) THEN 01608 DO xx=p_lb,coarse_gbo(1,1)-1 01609 x_att=coarse_gbo(1,1)+MODULO(xx-coarse_gbo(1,1),s(1)) 01610 IF (x_att>=my_coarse_bo(1,1).AND.x_att<=my_coarse_bo(2,1)) THEN 01611 send_size(p_old)=send_size(p_old)+coarse_slice_size 01612 END IF 01613 END DO 01614 DO xx=coarse_gbo(2,1)+1,p_ub 01615 x_att=coarse_gbo(1,1)+MODULO(xx-coarse_gbo(1,1),s(1)) 01616 IF (x_att>=my_coarse_bo(1,1).AND.x_att<=my_coarse_bo(2,1)) THEN 01617 send_size(p_old)=send_size(p_old)+coarse_slice_size 01618 END IF 01619 END DO 01620 END IF 01621 01622 p_old=p 01623 p_lb=FLOOR((x-2-f_shift(1))/2.) 01624 END IF 01625 END DO 01626 p_ub=FLOOR((fi_ub+3-f_shift(1))/2.) 01627 01628 send_size(p_old)=send_size(p_old)+(MIN(p_ub,my_coarse_bo(2,1)) & 01629 -MAX(p_lb,my_coarse_bo(1,1))+1)*coarse_slice_size 01630 01631 IF (pbc) THEN 01632 DO xx=p_lb,coarse_gbo(1,1)-1 01633 x_att=coarse_gbo(1,1)+MODULO(xx-coarse_gbo(1,1),s(1)) 01634 IF (x_att>=my_coarse_bo(1,1).AND.x_att<=my_coarse_bo(2,1)) THEN 01635 send_size(p_old)=send_size(p_old)+coarse_slice_size 01636 END IF 01637 END DO 01638 DO xx=coarse_gbo(2,1)+1,p_ub 01639 x_att=coarse_gbo(1,1)+MODULO(xx-coarse_gbo(1,1),s(1)) 01640 IF (x_att>=my_coarse_bo(1,1).AND.x_att<=my_coarse_bo(2,1)) THEN 01641 send_size(p_old)=send_size(p_old)+coarse_slice_size 01642 END IF 01643 END DO 01644 END IF 01645 ! ** offsets & alloc send-rcv 01646 01647 send_tot_size=0 01648 DO ip=0,n_procs-1 01649 send_offset(ip)=send_tot_size 01650 send_tot_size=send_tot_size+send_size(ip) 01651 END DO 01652 ALLOCATE(send_buf(0:send_tot_size-1),stat=stat) 01653 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 01654 01655 rcv_tot_size=0 01656 DO ip=0,n_procs-1 01657 rcv_offset(ip)=rcv_tot_size 01658 rcv_tot_size=rcv_tot_size+rcv_size(ip) 01659 END DO 01660 IF (.NOT.rcv_tot_size==(coarse_bo(2,1)-coarse_bo(1,1)+1)*coarse_slice_size) THEN 01661 CALL cp_assert(.FALSE.,& 01662 cp_failure_level,cp_assertion_failed,routineP,& 01663 "Error calculating rcv_tot_size "//& 01664 CPSourceFileRef,& 01665 error,failure) 01666 END IF 01667 ALLOCATE(rcv_buf(0:rcv_tot_size-1),stat=stat) 01668 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 01669 01670 ! ** fill send buffer 01671 01672 p_old=pos_of_x(fine_gbo(1,1)+MODULO(fi_lb-fine_gbo(1,1),sf)) 01673 p_lb=FLOOR((fi_lb-2-f_shift(1))/2.) 01674 sent_size=send_offset 01675 ss=my_coarse_bo(2,1)-my_coarse_bo(1,1)+1 01676 DO x=fi_lb,fi_ub 01677 p=pos_of_x(fine_gbo(1,1)+MODULO(x-fine_gbo(1,1),sf)) 01678 IF (p/=p_old) THEN 01679 shift=FLOOR((fine_gbo(1,1)+MODULO(x-1-fine_gbo(1,1),sf)-f_shift(1))/2._dp)-& 01680 FLOOR((x-1-f_shift(1))/2._dp) 01681 p_ub=FLOOR((x-1+3-f_shift(1))/2._dp) 01682 01683 IF (pbc) THEN 01684 DO xx=p_lb+shift,coarse_gbo(1,1)-1 01685 x_att=coarse_gbo(1,1)+MODULO(xx-coarse_gbo(1,1),sf) 01686 IF (x_att>=my_coarse_bo(1,1).AND.x_att<=my_coarse_bo(2,1)) THEN 01687 CALL dcopy(coarse_slice_size,& 01688 coarse_coeffs(x_att,my_coarse_bo(1,2),& 01689 my_coarse_bo(1,3)),ss,send_buf(sent_size(p_old)),1) 01690 sent_size(p_old)=sent_size(p_old)+coarse_slice_size 01691 END IF 01692 END DO 01693 END IF 01694 01695 ii=sent_size(p_old) 01696 DO k=coarse_bo(1,3),coarse_bo(2,3) 01697 DO j=coarse_bo(1,2),coarse_bo(2,2) 01698 DO i=MAX(p_lb+shift,my_coarse_bo(1,1)),MIN(p_ub+shift,my_coarse_bo(2,1)) 01699 send_buf(ii)=coarse_coeffs(i,j,k) 01700 ii=ii+1 01701 END DO 01702 END DO 01703 END DO 01704 sent_size(p_old)=ii 01705 01706 IF (pbc) THEN 01707 DO xx=coarse_gbo(2,1)+1,p_ub+shift 01708 x_att=coarse_gbo(1,1)+MODULO(xx-coarse_gbo(1,1),s(1)) 01709 IF (x_att>=my_coarse_bo(1,1).AND.x_att<=my_coarse_bo(2,1)) THEN 01710 CALL dcopy(coarse_slice_size,& 01711 coarse_coeffs(x_att,my_coarse_bo(1,2),& 01712 my_coarse_bo(1,3)),ss,& 01713 send_buf(sent_size(p_old)),1) 01714 sent_size(p_old)=sent_size(p_old)+coarse_slice_size 01715 END IF 01716 END DO 01717 END IF 01718 01719 p_old=p 01720 p_lb=FLOOR((x-2-f_shift(1))/2.) 01721 END IF 01722 END DO 01723 shift=FLOOR((fine_gbo(1,1)+MODULO(x-1-fine_gbo(1,1),sf)-f_shift(1))/2._dp)-& 01724 FLOOR((x-1-f_shift(1))/2._dp) 01725 p_ub=FLOOR((fi_ub+3-f_shift(1))/2.) 01726 01727 IF (pbc) THEN 01728 DO xx=p_lb+shift,coarse_gbo(1,1)-1 01729 x_att=coarse_gbo(1,1)+MODULO(xx-coarse_gbo(1,1),s(1)) 01730 IF (x_att>=my_coarse_bo(1,1).AND.x_att<=my_coarse_bo(2,1)) THEN 01731 CALL dcopy(coarse_slice_size,& 01732 coarse_coeffs(x_att,my_coarse_bo(1,2),& 01733 my_coarse_bo(1,3)),ss,send_buf(sent_size(p_old)),1) 01734 sent_size(p_old)=sent_size(p_old)+coarse_slice_size 01735 END IF 01736 END DO 01737 END IF 01738 01739 ii=sent_size(p_old) 01740 DO k=coarse_bo(1,3),coarse_bo(2,3) 01741 DO j=coarse_bo(1,2),coarse_bo(2,2) 01742 DO i=MAX(p_lb+shift,my_coarse_bo(1,1)),MIN(p_ub+shift,my_coarse_bo(2,1)) 01743 send_buf(ii)=coarse_coeffs(i,j,k) 01744 ii=ii+1 01745 END DO 01746 END DO 01747 END DO 01748 sent_size(p_old)=ii 01749 01750 IF (pbc) THEN 01751 DO xx=coarse_gbo(2,1)+1,p_ub+shift 01752 x_att=coarse_gbo(1,1)+MODULO(xx-coarse_gbo(1,1),s(1)) 01753 IF (x_att>=my_coarse_bo(1,1).AND.x_att<=my_coarse_bo(2,1)) THEN 01754 CALL dcopy(coarse_slice_size,& 01755 coarse_coeffs(x_att,my_coarse_bo(1,2),& 01756 my_coarse_bo(1,3)),ss,send_buf(sent_size(p_old)),1) 01757 sent_size(p_old)=sent_size(p_old)+coarse_slice_size 01758 END IF 01759 END DO 01760 END IF 01761 01762 CPPostcondition(ALL(sent_size(:n_procs-2)==send_offset(1:)),cp_failure_level,routineP,error,failure) 01763 CPPostcondition(sent_size(n_procs-1)==send_tot_size,cp_failure_level,routineP,error,failure) 01764 ! test send/rcv sizes 01765 CALL mp_alltoall(send_size,real_rcv_size,1,coarse_coeffs_pw%pw_grid%para%group) 01766 CPAssert(ALL(real_rcv_size==rcv_size),cp_failure_level,routineP,error,failure) 01767 ! all2all 01768 CALL mp_alltoall( sb=send_buf, scount=send_size, sdispl=send_offset,& 01769 rb=rcv_buf, rcount=rcv_size, rdispl=rcv_offset, & 01770 group=coarse_coeffs_pw%pw_grid%para%group ) 01771 01772 ! ** reorder rcv buffer 01773 ! (actually reordering should be needed only with pbc) 01774 01775 ALLOCATE(coarse_coeffs(coarse_bo(1,1):coarse_bo(2,1),& 01776 coarse_bo(1,2):coarse_bo(2,2),& 01777 coarse_bo(1,3):coarse_bo(2,3)),stat=stat) 01778 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 01779 01780 my_lb=MAX(coarse_gbo(1,1),coarse_bo(1,1)) 01781 my_ub=MIN(coarse_gbo(2,1),coarse_bo(2,1)) 01782 pos_of_x => coarse_coeffs_pw%pw_grid%para%pos_of_x 01783 sent_size=rcv_offset 01784 ss=coarse_bo(2,1)-coarse_bo(1,1)+1 01785 DO x=my_ub+1,coarse_bo(2,1) 01786 p_old=pos_of_x(coarse_gbo(1,1)+MODULO(x-coarse_gbo(1,1),s(1))) 01787 CALL dcopy(coarse_slice_size,& 01788 rcv_buf(sent_size(p_old)),1,& 01789 coarse_coeffs(x,coarse_bo(1,2),& 01790 coarse_bo(1,3)),ss) 01791 sent_size(p_old)=sent_size(p_old)+coarse_slice_size 01792 END DO 01793 p_old=pos_of_x(coarse_gbo(1,1) & 01794 +MODULO(my_lb-coarse_gbo(1,1),s(1))) 01795 p_lb=my_lb 01796 DO x=my_lb,my_ub 01797 p=pos_of_x(coarse_gbo(1,1)+MODULO(x-coarse_gbo(1,1),s(1))) 01798 IF (p/=p_old) THEN 01799 p_ub=x-1 01800 01801 ii=sent_size(p_old) 01802 DO k=coarse_bo(1,3),coarse_bo(2,3) 01803 DO j=coarse_bo(1,2),coarse_bo(2,2) 01804 DO i=p_lb,p_ub 01805 coarse_coeffs(i,j,k)=rcv_buf(ii) 01806 ii=ii+1 01807 END DO 01808 END DO 01809 END DO 01810 sent_size(p_old)=ii 01811 01812 p_lb=x 01813 p_old=p 01814 END IF 01815 rcv_size(p)=rcv_size(p)+coarse_slice_size 01816 END DO 01817 p_ub=my_ub 01818 ii=sent_size(p_old) 01819 DO k=coarse_bo(1,3),coarse_bo(2,3) 01820 DO j=coarse_bo(1,2),coarse_bo(2,2) 01821 DO i=p_lb,p_ub 01822 coarse_coeffs(i,j,k)=rcv_buf(ii) 01823 ii=ii+1 01824 END DO 01825 END DO 01826 END DO 01827 sent_size(p_old)=ii 01828 DO x=coarse_bo(1,1),my_lb-1 01829 p_old=pos_of_x(coarse_gbo(1,1)+MODULO(x-coarse_gbo(1,1),s(1))) 01830 CALL dcopy(coarse_slice_size,& 01831 rcv_buf(sent_size(p_old)),1,& 01832 coarse_coeffs(x,coarse_bo(1,2),& 01833 coarse_bo(1,3)),ss) 01834 sent_size(p_old)=sent_size(p_old)+coarse_slice_size 01835 END DO 01836 01837 CPPostcondition(ALL(sent_size(0:n_procs-2)==rcv_offset(1:)),cp_failure_level,routineP,error,failure) 01838 CPPostcondition(sent_size(n_procs-1)==rcv_tot_size,cp_failure_level,routineP,error,failure) 01839 01840 ! dealloc 01841 DEALLOCATE(send_size,send_offset, rcv_size, rcv_offset, stat=stat) 01842 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 01843 DEALLOCATE(send_buf,rcv_buf,real_rcv_size,stat=stat) 01844 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 01845 CALL timestop(handle2) 01846 01847 END IF 01848 fine_values => fine_values_pw%cr3d 01849 w_0=(/ weights_1d(3),weights_1d(1),weights_1d(3) /) 01850 w_1=(/ weights_1d(4),weights_1d(2),weights_1d(2),weights_1d(4) /) 01851 01852 DO k=coarse_bo(1,3),coarse_bo(2,3) 01853 DO ik=-3,3 01854 IF (pbc) THEN 01855 wk=weights_1d(ABS(ik)+1) 01856 fk=fine_gbo(1,3)+MODULO(2*k+ik-fine_gbo(1,3)+f_shift(3),2*s(3)) 01857 ELSE 01858 fk=2*k+ik+f_shift(3) 01859 IF (fk<=fine_bo(1,3)+1.OR.fk>=fine_bo(2,3)-1) THEN 01860 IF (fk<fine_bo(1,3).OR.fk>fine_bo(2,3)) CYCLE 01861 IF (fk==fine_bo(1,3).OR.fk==fine_bo(2,3)) THEN 01862 IF (ik/=0) CYCLE 01863 wk=w_border0 01864 ELSE IF (fk==2*coarse_bo(1,3)+1+f_shift(3)) THEN 01865 SELECT CASE(ik) 01866 CASE(1) 01867 wk=w_border1(1) 01868 CASE(-1) 01869 wk=w_border1(2) 01870 CASE(-3) 01871 wk=w_border1(3) 01872 CASE default 01873 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 01874 CYCLE 01875 END SELECT 01876 ELSE 01877 SELECT CASE(ik) 01878 CASE(3) 01879 wk=w_border1(3) 01880 CASE(1) 01881 wk=w_border1(2) 01882 CASE(-1) 01883 wk=w_border1(1) 01884 CASE default 01885 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 01886 CYCLE 01887 END SELECT 01888 END IF 01889 ELSE 01890 wk=weights_1d(ABS(ik)+1) 01891 END IF 01892 END IF 01893 DO j=coarse_bo(1,2),coarse_bo(2,2) 01894 DO ij=-3,3 01895 IF (pbc) THEN 01896 wj=weights_1d(ABS(ij)+1)*wk 01897 fj=fine_gbo(1,2)+MODULO(2*j+ij-fine_gbo(1,2)+f_shift(2),2*s(2)) 01898 ELSE 01899 fj=2*j+ij+f_shift(2) 01900 IF (fj<=fine_bo(1,2)+1.OR.fj>=fine_bo(2,2)-1) THEN 01901 IF (fj<fine_bo(1,2).OR.fj>fine_bo(2,2)) CYCLE 01902 IF (fj==fine_bo(1,2).OR.fj==fine_bo(2,2)) THEN 01903 IF (ij/=0) CYCLE 01904 wj=w_border0*wk 01905 ELSE IF (fj==2*coarse_bo(1,2)+1+f_shift(2)) THEN 01906 SELECT CASE(ij) 01907 CASE(1) 01908 wj=w_border1(1)*wk 01909 CASE(-1) 01910 wj=w_border1(2)*wk 01911 CASE(-3) 01912 wj=w_border1(3)*wk 01913 CASE default 01914 CYCLE 01915 END SELECT 01916 ELSE 01917 SELECT CASE(ij) 01918 CASE(-1) 01919 wj=w_border1(1)*wk 01920 CASE(1) 01921 wj=w_border1(2)*wk 01922 CASE(3) 01923 wj=w_border1(3)*wk 01924 CASE default 01925 CYCLE 01926 END SELECT 01927 END IF 01928 ELSE 01929 wj=weights_1d(ABS(ij)+1)*wk 01930 END IF 01931 END IF 01932 01933 IF (fine_bo(2,1)-fine_bo(1,1)<7.or.safe_calc) THEN 01934 ! CALL timeset(routineN//"_safe",handle2) 01935 DO i=coarse_bo(1,1),coarse_bo(2,1) 01936 DO ii=-3,3 01937 IF (pbc.and..not.is_split) THEN 01938 wi=weights_1d(ABS(ii)+1)*wj 01939 fi=fine_gbo(1,1)+MODULO(2*i+ii-fine_gbo(1,1)+f_shift(1),2*s(1)) 01940 ELSE 01941 fi=2*i+ii+f_shift(1) 01942 IF (fi<fine_bo(1,1).OR.fi>fine_bo(2,1)) CYCLE 01943 IF (.NOT.pbc.AND.(fi<=fine_gbo(1,1)+1.OR. & 01944 fi>=fine_gbo(2,1)-1)) THEN 01945 IF (fi==fine_gbo(1,1).OR.fi==fine_gbo(2,1)) THEN 01946 IF (ii/=0) CYCLE 01947 wi=w_border0*wj 01948 ELSE IF (fi==fine_gbo(1,1)+1) THEN 01949 SELECT CASE(ii) 01950 CASE(1) 01951 wi=w_border1(1)*wj 01952 CASE(-1) 01953 wi=w_border1(2)*wj 01954 CASE(-3) 01955 wi=w_border1(3)*wj 01956 CASE default 01957 CYCLE 01958 END SELECT 01959 ELSE 01960 SELECT CASE(ii) 01961 CASE(-1) 01962 wi=w_border1(1)*wj 01963 CASE(1) 01964 wi=w_border1(2)*wj 01965 CASE(3) 01966 wi=w_border1(3)*wj 01967 CASE default 01968 CYCLE 01969 END SELECT 01970 END IF 01971 ELSE 01972 wi=weights_1d(ABS(ii)+1)*wj 01973 END IF 01974 END IF 01975 fine_values(fi,fj,fk)=& 01976 fine_values(fi,fj,fk)+& 01977 wi*coarse_coeffs(i,j,k) 01978 END DO 01979 END DO 01980 ! CALL timestop(handle2) 01981 ELSE 01982 ! CALL timeset(routineN//"_core1",handle2) 01983 ww0=wj*w_0 01984 ww1=wj*w_1 01985 IF (pbc.AND..NOT.is_split) THEN 01986 v3=coarse_coeffs(coarse_bo(2,1),j,k) 01987 i=coarse_bo(1,1) 01988 fi=2*i+f_shift(1) 01989 v0= coarse_coeffs(i,j,k) 01990 v1= coarse_coeffs(i+1,j,k) 01991 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 01992 ww0(1)*v3+ww0(2)*v0+ww0(3)*v1 01993 v2= coarse_coeffs(i+2,j,k) 01994 fi=fi+1 01995 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 01996 ww1(1)*v3+ww1(2)*v0+ww1(3)*v1+ww1(4)*v2 01997 ELSE IF (.NOT.has_i_lbound) THEN 01998 i=coarse_bo(1,1) 01999 fi=2*i+f_shift(1) 02000 v0= coarse_coeffs(i,j,k) 02001 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02002 w_border0*wj*v0 02003 v1= coarse_coeffs(i+1,j,k) 02004 v2= coarse_coeffs(i+2,j,k) 02005 fi=fi+1 02006 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02007 wj*(w_border1(1)*v0+w_border1(2)*v1+& 02008 w_border1(3)*v2) 02009 ELSE 02010 i=coarse_bo(1,1) 02011 v0= coarse_coeffs(i,j,k) 02012 v1= coarse_coeffs(i+1,j,k) 02013 v2= coarse_coeffs(i+2,j,k) 02014 fi=2*i+f_shift(1)+1 02015 IF (.NOT.(fi+1==fine_bo(1,1).OR.& 02016 fi+2==fine_bo(1,1))) THEN 02017 CALL cp_assert(.FALSE.,cp_failure_level,& 02018 cp_assertion_failed,routineN,& 02019 "unexpected start index "//& 02020 TRIM(cp_to_string(coarse_bo(1,1)))//" "//& 02021 TRIM(cp_to_string(fi)),error,failure) 02022 END IF 02023 END IF 02024 fi=fi+1 02025 IF (fi>=fine_bo(1,1)) THEN 02026 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02027 ww0(1)*v0+ww0(2)*v1+& 02028 ww0(3)*v2 02029 ELSE 02030 CPPrecondition(fi+1==fine_bo(1,1),cp_failure_level,routineP,error,failure) 02031 END IF 02032 ! CALL timestop(handle2) 02033 ! CALL timeset(routineN//"_core",handle2) 02034 DO i=coarse_bo(1,1)+3,FLOOR((fine_bo(2,1)-f_shift(1))/2.)-3,4 02035 v3=coarse_coeffs(i,j,k) 02036 fi=fi+1 02037 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02038 (ww1(1)*v0+ww1(2)*v1+& 02039 ww1(3)*v2+ww1(4)*v3) 02040 fi=fi+1 02041 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02042 (ww0(1)*v1+ww0(2)*v2+& 02043 ww0(3)*v3) 02044 v0=coarse_coeffs(i+1,j,k) 02045 fi=fi+1 02046 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02047 (ww1(4)*v0+ww1(1)*v1+& 02048 ww1(2)*v2+ww1(3)*v3) 02049 fi=fi+1 02050 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02051 (ww0(1)*v2+ww0(2)*v3+& 02052 ww0(3)*v0) 02053 v1=coarse_coeffs(i+2,j,k) 02054 fi=fi+1 02055 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02056 (ww1(3)*v0+ww1(4)*v1+& 02057 ww1(1)*v2+ww1(2)*v3) 02058 fi=fi+1 02059 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02060 (ww0(1)*v3+ww0(2)*v0+& 02061 ww0(3)*v1) 02062 v2=coarse_coeffs(i+3,j,k) 02063 fi=fi+1 02064 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02065 (ww1(2)*v0+ww1(3)*v1+& 02066 ww1(4)*v2+ww1(1)*v3) 02067 fi=fi+1 02068 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02069 (ww0(1)*v0+ww0(2)*v1+& 02070 ww0(3)*v2) 02071 END DO 02072 ! CALL timestop(handle2) 02073 ! CALL timeset(routineN//"_clean",handle2) 02074 rest_b=MODULO(FLOOR((fine_bo(2,1)-f_shift(1))/2.)-coarse_bo(1,1)-3+1,4) 02075 IF (rest_b>0) THEN 02076 i=FLOOR((fine_bo(2,1)-f_shift(1))/2.)-rest_b+1 02077 v3=coarse_coeffs(i,j,k) 02078 fi=fi+1 02079 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02080 (ww1(1)*v0+ww1(2)*v1+& 02081 ww1(3)*v2+ww1(4)*v3) 02082 fi=fi+1 02083 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02084 (ww0(1)*v1+ww0(2)*v2+& 02085 ww0(3)*v3) 02086 IF (rest_b>1) THEN 02087 v0=coarse_coeffs(i+1,j,k) 02088 fi=fi+1 02089 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02090 (ww1(4)*v0+ww1(1)*v1+& 02091 ww1(2)*v2+ww1(3)*v3) 02092 fi=fi+1 02093 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02094 (ww0(1)*v2+ww0(2)*v3+& 02095 ww0(3)*v0) 02096 IF (rest_b>2) THEN 02097 v1=coarse_coeffs(i+2,j,k) 02098 fi=fi+1 02099 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02100 (ww1(3)*v0+ww1(4)*v1+& 02101 ww1(1)*v2+ww1(2)*v3) 02102 fi=fi+1 02103 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02104 (ww0(1)*v3+ww0(2)*v0+& 02105 ww0(3)*v1) 02106 IF (pbc.and..not.is_split) THEN 02107 v2=coarse_coeffs(coarse_bo(1,1),j,k) 02108 fi=fi+1 02109 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02110 ww1(1)*v3+ww1(2)*v0+ww1(3)*v1+ww1(4)*v2 02111 fi=fi+1 02112 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02113 ww0(1)*v0+ww0(2)*v1+ww0(3)*v2 02114 v3=coarse_coeffs(coarse_bo(1,1)+1,j,k) 02115 fi=fi+1 02116 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02117 ww1(1)*v0+ww1(2)*v1+ww1(3)*v2+ww1(4)*v3 02118 ELSE IF (has_i_ubound) THEN 02119 v2=coarse_coeffs(i+3,j,k) 02120 fi=fi+1 02121 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02122 ww1(1)*v3+ww1(2)*v0+ww1(3)*v1+ww1(4)*v2 02123 fi=fi+1 02124 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02125 ww0(1)*v0+ww0(2)*v1+ww0(3)*v2 02126 IF (fi+1==fine_bo(2,1)) THEN 02127 v3=coarse_coeffs(i+4,j,k) 02128 fi=fi+1 02129 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02130 ww1(1)*v0+ww1(2)*v1+ww1(3)*v2+ww1(4)*v3 02131 END IF 02132 ELSE 02133 fi=fi+1 02134 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02135 wj*(w_border1(3)*v3+w_border1(2)*v0+& 02136 w_border1(1)*v1) 02137 fi=fi+1 02138 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02139 w_border0*wj*v1 02140 END IF 02141 ELSE IF (pbc.and..not.is_split) THEN 02142 v1=coarse_coeffs(coarse_bo(1,1),j,k) 02143 fi=fi+1 02144 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02145 ww1(1)*v2+ww1(2)*v3+ww1(3)*v0+ww1(4)*v1 02146 fi=fi+1 02147 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02148 ww0(1)*v3+ww0(2)*v0+ww0(3)*v1 02149 v2=coarse_coeffs(coarse_bo(1,1)+1,j,k) 02150 fi=fi+1 02151 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02152 ww1(1)*v3+ww1(2)*v0+ww1(3)*v1+ww1(4)*v2 02153 ELSE IF(has_i_ubound) THEN 02154 v1=coarse_coeffs(i+2,j,k) 02155 fi=fi+1 02156 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02157 ww1(1)*v2+ww1(2)*v3+ww1(3)*v0+ww1(4)*v1 02158 fi=fi+1 02159 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02160 ww0(1)*v3+ww0(2)*v0+ww0(3)*v1 02161 IF (fi+1==fine_bo(2,1)) THEN 02162 v2=coarse_coeffs(i+3,j,k) 02163 fi=fi+1 02164 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02165 ww1(1)*v3+ww1(2)*v0+ww1(3)*v1+ww1(4)*v2 02166 END IF 02167 ELSE 02168 fi=fi+1 02169 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02170 wj*(w_border1(3)*v2+w_border1(2)*v3+& 02171 w_border1(1)*v0) 02172 fi=fi+1 02173 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02174 w_border0*wj*v0 02175 END IF 02176 ELSE IF (pbc.and..not.is_split) THEN 02177 v0=coarse_coeffs(coarse_bo(1,1),j,k) 02178 fi=fi+1 02179 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02180 ww1(1)*v1+ww1(2)*v2+ww1(3)*v3+ww1(4)*v0 02181 fi=fi+1 02182 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02183 ww0(1)*v2+ww0(2)*v3+ww0(3)*v0 02184 v1=coarse_coeffs(coarse_bo(1,1)+1,j,k) 02185 fi=fi+1 02186 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02187 ww1(1)*v2+ww1(2)*v3+ww1(3)*v0+ww1(4)*v1 02188 ELSE IF (has_i_ubound) THEN 02189 v0=coarse_coeffs(i+1,j,k) 02190 fi=fi+1 02191 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02192 ww1(1)*v1+ww1(2)*v2+ww1(3)*v3+ww1(4)*v0 02193 fi=fi+1 02194 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02195 ww0(1)*v2+ww0(2)*v3+ww0(3)*v0 02196 IF (fi+1==fine_bo(2,1)) THEN 02197 v1=coarse_coeffs(i+2,j,k) 02198 fi=fi+1 02199 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02200 ww1(1)*v2+ww1(2)*v3+ww1(3)*v0+ww1(4)*v1 02201 END IF 02202 ELSE 02203 fi=fi+1 02204 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02205 wj*(w_border1(3)*v1+w_border1(2)*v2+& 02206 w_border1(1)*v3) 02207 fi=fi+1 02208 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02209 w_border0*wj*v3 02210 END IF 02211 ELSE IF (pbc.and..not.is_split) THEN 02212 v3=coarse_coeffs(coarse_bo(1,1),j,k) 02213 fi=fi+1 02214 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02215 ww1(1)*v0+ww1(2)*v1+ww1(3)*v2+ww1(4)*v3 02216 fi=fi+1 02217 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02218 ww0(1)*v1+ww0(2)*v2+ww0(3)*v3 02219 v0=coarse_coeffs(coarse_bo(1,1)+1,j,k) 02220 fi=fi+1 02221 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02222 ww1(1)*v1+ww1(2)*v2+ww1(3)*v3+ww1(4)*v0 02223 ELSE IF (has_i_ubound) THEN 02224 v3=coarse_coeffs(i,j,k) 02225 fi=fi+1 02226 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02227 ww1(1)*v0+ww1(2)*v1+ww1(3)*v2+ww1(4)*v3 02228 fi=fi+1 02229 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02230 ww0(1)*v1+ww0(2)*v2+ww0(3)*v3 02231 IF (fi+1==fine_bo(2,1)) THEN 02232 v0=coarse_coeffs(i+1,j,k) 02233 fi=fi+1 02234 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02235 ww1(1)*v1+ww1(2)*v2+ww1(3)*v3+ww1(4)*v0 02236 END IF 02237 ELSE 02238 fi=fi+1 02239 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02240 wj*(w_border1(3)*v0+w_border1(2)*v1+& 02241 w_border1(1)*v2) 02242 fi=fi+1 02243 fine_values(fi,fj,fk)=fine_values(fi,fj,fk)+& 02244 w_border0*wj*v2 02245 END IF 02246 CPPostcondition(fi==fine_bo(2,1),cp_failure_level,routineP,error,failure) 02247 END IF 02248 ! CALL timestop(handle2) 02249 END DO 02250 END DO 02251 END DO 02252 END DO 02253 02254 IF (is_split) THEN 02255 DEALLOCATE(coarse_coeffs,stat=stat) 02256 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02257 END IF 02258 END IF 02259 CALL timestop(handle) 02260 END SUBROUTINE add_coarse2fine 02261 02262 ! ***************************************************************************** 02299 SUBROUTINE add_fine2coarse(fine_values_pw,coarse_coeffs_pw,& 02300 weights_1d,w_border0,w_border1,pbc,safe_computation,error) 02301 TYPE(pw_type), POINTER :: fine_values_pw, 02302 coarse_coeffs_pw 02303 REAL(kind=dp), DIMENSION(4), INTENT(in) :: weights_1d 02304 REAL(kind=dp), INTENT(in) :: w_border0 02305 REAL(kind=dp), DIMENSION(3), INTENT(in) :: w_border1 02306 LOGICAL, INTENT(in) :: pbc 02307 LOGICAL, INTENT(in), OPTIONAL :: safe_computation 02308 TYPE(cp_error_type), INTENT(inout) :: error 02309 02310 CHARACTER(len=*), PARAMETER :: routineN = 'add_fine2coarse', 02311 routineP = moduleN//':'//routineN 02312 02313 INTEGER :: coarse_slice_size, f_shift(3), fi, fj, fk, handle, handle2, i, 02314 ii, ij, ik, ip, j, k, n_procs, p, p_old, rcv_tot_size, rest_b, s(3), 02315 send_tot_size, ss, stat, x, x_att 02316 INTEGER, ALLOCATABLE, DIMENSION(:) :: pp_lb, pp_ub, rcv_offset, 02317 rcv_size, real_rcv_size, 02318 send_offset, send_size, 02319 sent_size 02320 INTEGER, DIMENSION(2, 3) :: coarse_bo, coarse_gbo, 02321 fine_bo, fine_gbo, 02322 my_coarse_bo 02323 INTEGER, DIMENSION(:), POINTER :: pos_of_x 02324 LOGICAL :: failure, has_i_lbound, 02325 has_i_ubound, is_split, 02326 local_data, safe_calc 02327 REAL(kind=dp) :: vv0, vv1, vv2, vv3, vv4, vv5, 02328 vv6, vv7, wi, wj, wk 02329 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: rcv_buf, send_buf 02330 REAL(kind=dp), DIMENSION(3) :: w_0, ww0 02331 REAL(kind=dp), DIMENSION(4) :: w_1, ww1 02332 REAL(kind=dp), DIMENSION(:, :, :), 02333 POINTER :: coarse_coeffs, fine_values 02334 02335 failure=.FALSE. 02336 CALL timeset(routineN,handle) 02337 02338 IF (.NOT. failure) THEN 02339 safe_calc=.FALSE. 02340 IF (PRESENT(safe_computation)) safe_calc=safe_computation 02341 02342 my_coarse_bo=coarse_coeffs_pw%pw_grid%bounds_local 02343 coarse_gbo=coarse_coeffs_pw%pw_grid%bounds 02344 fine_bo=fine_values_pw%pw_grid%bounds_local 02345 fine_gbo=fine_values_pw%pw_grid%bounds 02346 f_shift=fine_gbo(1,:)-2*coarse_gbo(1,:) 02347 is_split= ANY(coarse_gbo(:,1) /= my_coarse_bo(:,1)) 02348 coarse_bo=my_coarse_bo 02349 IF (fine_bo(1,1)<=fine_bo(2,1)) THEN 02350 coarse_bo(1,1)=FLOOR(REAL(fine_bo(1,1)-f_shift(1),dp)/2._dp)-1 02351 coarse_bo(2,1)=FLOOR(REAL(fine_bo(2,1)+1-f_shift(1),dp)/2._dp)+1 02352 ELSE 02353 coarse_bo(1,1)=coarse_gbo(2,1) 02354 coarse_bo(2,1)=coarse_gbo(2,1)-1 02355 END IF 02356 IF (.NOT.is_split.OR..NOT.pbc) THEN 02357 coarse_bo(1,1)=MAX(coarse_gbo(1,1),coarse_bo(1,1)) 02358 coarse_bo(2,1)=MIN(coarse_gbo(2,1),coarse_bo(2,1)) 02359 END IF 02360 has_i_ubound=(fine_gbo(2,1)/=fine_bo(2,1)).or.pbc.and.is_split 02361 has_i_lbound=(fine_gbo(1,1)/=fine_bo(1,1)).or.pbc.and.is_split 02362 02363 IF (pbc) THEN 02364 CPPrecondition(ALL(fine_gbo(1,:)==2*coarse_gbo(1,:)+f_shift),cp_failure_level,routineP,error,failure) 02365 CPPrecondition(ALL(fine_gbo(2,:)==2*coarse_gbo(2,:)+f_shift+1),cp_failure_level,routineP,error,failure) 02366 ELSE 02367 CPPrecondition(ALL(fine_gbo(2,:)==2*coarse_gbo(2,:)+f_shift),cp_failure_level,routineP,error,failure) 02368 CPPrecondition(ALL(fine_gbo(1,:)==2*coarse_gbo(1,:)+f_shift),cp_failure_level,routineP,error,failure) 02369 END IF 02370 CPPrecondition(coarse_gbo(2,1)-coarse_gbo(1,2)>1,cp_failure_level,routineP,error,failure) 02371 local_data=is_split ! ANY(coarse_bo/=my_coarse_bo) 02372 IF (local_data) THEN 02373 ALLOCATE(coarse_coeffs(coarse_bo(1,1):coarse_bo(2,1),& 02374 coarse_bo(1,2):coarse_bo(2,2),& 02375 coarse_bo(1,3):coarse_bo(2,3)),stat=stat) 02376 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02377 coarse_coeffs=0._dp 02378 ELSE 02379 coarse_coeffs => coarse_coeffs_pw%cr3d 02380 END IF 02381 02382 fine_values => fine_values_pw%cr3d 02383 w_0=(/ weights_1d(3),weights_1d(1),weights_1d(3) /) 02384 w_1=(/ weights_1d(4),weights_1d(2),weights_1d(2),weights_1d(4) /) 02385 02386 DO i=1,3 02387 s(i)=coarse_gbo(2,i)-coarse_gbo(1,i)+1 02388 END DO 02389 IF (ANY(s<1)) RETURN 02390 02391 DO k=coarse_bo(1,3),coarse_bo(2,3) 02392 DO ik=-3,3 02393 IF (pbc) THEN 02394 wk=weights_1d(ABS(ik)+1) 02395 fk=fine_gbo(1,3)+MODULO(2*k+ik-fine_gbo(1,3)+f_shift(3),2*s(3)) 02396 ELSE 02397 fk=2*k+ik+f_shift(3) 02398 IF (fk<=fine_bo(1,3)+1.OR.fk>=fine_bo(2,3)-1) THEN 02399 IF (fk<fine_bo(1,3).OR.fk>fine_bo(2,3)) CYCLE 02400 IF (fk==fine_bo(1,3).OR.fk==fine_bo(2,3)) THEN 02401 IF (ik/=0) CYCLE 02402 wk=w_border0 02403 ELSE IF (fk==fine_bo(1,3)+1) THEN 02404 SELECT CASE(ik) 02405 CASE(1) 02406 wk=w_border1(1) 02407 CASE(-1) 02408 wk=w_border1(2) 02409 CASE(-3) 02410 wk=w_border1(3) 02411 CASE default 02412 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 02413 CYCLE 02414 END SELECT 02415 ELSE 02416 SELECT CASE(ik) 02417 CASE(3) 02418 wk=w_border1(3) 02419 CASE(1) 02420 wk=w_border1(2) 02421 CASE(-1) 02422 wk=w_border1(1) 02423 CASE default 02424 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 02425 CYCLE 02426 END SELECT 02427 END IF 02428 ELSE 02429 wk=weights_1d(ABS(ik)+1) 02430 END IF 02431 END IF 02432 DO j=coarse_bo(1,2),coarse_bo(2,2) 02433 DO ij=-3,3 02434 IF (pbc) THEN 02435 fj=fine_gbo(1,2)+MODULO(2*j+ij-fine_gbo(1,2)+f_shift(2),& 02436 2*s(2)) 02437 wj=weights_1d(ABS(ij)+1)*wk 02438 ELSE 02439 fj=2*j+ij+f_shift(2) 02440 IF (fj<=fine_bo(1,2)+1.OR.fj>=fine_bo(2,2)-1) THEN 02441 IF (fj<fine_bo(1,2).OR.fj>fine_bo(2,2)) CYCLE 02442 IF (fj==fine_bo(1,2).OR.fj==fine_bo(2,2)) THEN 02443 IF (ij/=0) CYCLE 02444 wj=w_border0*wk 02445 ELSE IF (fj==fine_bo(1,2)+1) THEN 02446 SELECT CASE(ij) 02447 CASE(1) 02448 wj=w_border1(1)*wk 02449 CASE(-1) 02450 wj=w_border1(2)*wk 02451 CASE(-3) 02452 wj=w_border1(3)*wk 02453 CASE default 02454 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 02455 CYCLE 02456 END SELECT 02457 ELSE 02458 SELECT CASE(ij) 02459 CASE(-1) 02460 wj=w_border1(1)*wk 02461 CASE(1) 02462 wj=w_border1(2)*wk 02463 CASE(3) 02464 wj=w_border1(3)*wk 02465 CASE default 02466 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 02467 CYCLE 02468 END SELECT 02469 END IF 02470 ELSE 02471 wj=weights_1d(ABS(ij)+1)*wk 02472 END IF 02473 END IF 02474 02475 IF (coarse_bo(2,1)-coarse_bo(1,1)<7.or.safe_calc) THEN 02476 DO i=coarse_bo(1,1),coarse_bo(2,1) 02477 DO ii=-3,3 02478 IF (pbc.AND..NOT.is_split) THEN 02479 wi=weights_1d(ABS(ii)+1)*wj 02480 fi=fine_gbo(1,1)+MODULO(2*i+ii-fine_gbo(1,1)+f_shift(1),2*s(1)) 02481 ELSE 02482 fi=2*i+ii+f_shift(1) 02483 IF (fi<fine_bo(1,1).OR.fi>fine_bo(2,1)) CYCLE 02484 IF (((.NOT.pbc).AND.fi<=fine_gbo(1,1)+1).OR. & 02485 ((.NOT.pbc).AND.fi>=fine_gbo(2,1)-1)) THEN 02486 IF (fi==fine_gbo(1,1).OR.fi==fine_gbo(2,1)) THEN 02487 IF (ii/=0) CYCLE 02488 wi=w_border0*wj 02489 ELSE IF (fi==fine_gbo(1,1)+1) THEN 02490 SELECT CASE(ii) 02491 CASE(1) 02492 wi=w_border1(1)*wj 02493 CASE(-1) 02494 wi=w_border1(2)*wj 02495 CASE(-3) 02496 wi=w_border1(3)*wj 02497 CASE default 02498 CYCLE 02499 END SELECT 02500 ELSE 02501 SELECT CASE(ii) 02502 CASE(-1) 02503 wi=w_border1(1)*wj 02504 CASE(1) 02505 wi=w_border1(2)*wj 02506 CASE(3) 02507 wi=w_border1(3)*wj 02508 CASE default 02509 CYCLE 02510 END SELECT 02511 END IF 02512 ELSE 02513 wi=weights_1d(ABS(ii)+1)*wj 02514 END IF 02515 END IF 02516 coarse_coeffs(i,j,k)=& 02517 coarse_coeffs(i,j,k)+& 02518 wi*fine_values(fi,fj,fk) 02519 END DO 02520 END DO 02521 ELSE 02522 ww0=wj*w_0 02523 ww1=wj*w_1 02524 IF (pbc.AND..NOT.is_split) THEN 02525 i=coarse_bo(1,1)-1 02526 vv2=fine_values(fine_bo(2,1)-2,fj,fk) 02527 vv3=fine_values(fine_bo(2,1)-1,fj,fk) 02528 vv4=fine_values(fine_bo(2,1),fj,fk) 02529 fi=fine_bo(1,1) 02530 vv5=fine_values(fi,fj,fk) 02531 fi=fi+1 02532 vv6=fine_values(fi,fj,fk) 02533 fi=fi+1 02534 vv7=fine_values(fi,fj,fk) 02535 coarse_coeffs(i+1,j,k)=coarse_coeffs(i+1,j,k)& 02536 +ww1(4)*vv2+ww0(3)*vv3+ww1(3)*vv4+ww0(2)*vv5+ww1(2)*vv6+ww0(1)*vv7 02537 coarse_coeffs(i+2,j,k)=coarse_coeffs(i+2,j,k)& 02538 +ww1(4)*vv4+ww0(3)*vv5+ww1(3)*vv6+ww0(2)*vv7 02539 coarse_coeffs(i+3,j,k)=coarse_coeffs(i+3,j,k)& 02540 +ww1(4)*vv6+ww0(3)*vv7 02541 ELSE IF (has_i_lbound) THEN 02542 i=coarse_bo(1,1) 02543 fi=fine_bo(1,1)-1 02544 IF (i+1==FLOOR((fine_bo(1,1)+1-f_shift(1))/2._dp)) THEN 02545 fi=fi+1 02546 vv0=fine_values(fi,fj,fk) 02547 coarse_coeffs(i,j,k)=coarse_coeffs(i,j,k)+& 02548 02549 vv0*ww0(3) 02550 coarse_coeffs(i+1,j,k)=coarse_coeffs(i+1,j,k)+& 02551 02552 vv0*ww0(2) 02553 coarse_coeffs(i+2,j,k)=coarse_coeffs(i+2,j,k)+& 02554 02555 vv0*ww0(1) 02556 END IF 02557 ELSE 02558 i=coarse_bo(1,1) 02559 fi=2*i+f_shift(1) 02560 vv0=fine_values(fi,fj,fk) 02561 fi=fi+1 02562 vv1=fine_values(fi,fj,fk) 02563 fi=fi+1 02564 vv2=fine_values(fi,fj,fk) 02565 coarse_coeffs(i,j,k)=coarse_coeffs(i,j,k)+& 02566 (vv0*w_border0+vv1*w_border1(1))*wj+vv2*ww0(1) 02567 coarse_coeffs(i+1,j,k)=coarse_coeffs(i+1,j,k)+& 02568 wj*w_border1(2)*vv1+ww0(2)*vv2 02569 coarse_coeffs(i+2,j,k)=coarse_coeffs(i+2,j,k)+& 02570 wj*w_border1(3)*vv1+ww0(3)*vv2 02571 END IF 02572 DO i=coarse_bo(1,1)+3,FLOOR((fine_bo(2,1)-f_shift(1))/2._dp)-3,4 02573 fi=fi+1 02574 vv0=fine_values(fi,fj,fk) 02575 fi=fi+1 02576 vv1=fine_values(fi,fj,fk) 02577 fi=fi+1 02578 vv2=fine_values(fi,fj,fk) 02579 fi=fi+1 02580 vv3=fine_values(fi,fj,fk) 02581 fi=fi+1 02582 vv4=fine_values(fi,fj,fk) 02583 fi=fi+1 02584 vv5=fine_values(fi,fj,fk) 02585 fi=fi+1 02586 vv6=fine_values(fi,fj,fk) 02587 fi=fi+1 02588 vv7=fine_values(fi,fj,fk) 02589 coarse_coeffs(i-3,j,k)=coarse_coeffs(i-3,j,k)& 02590 +ww1(1)*vv0 02591 coarse_coeffs(i-2,j,k)=coarse_coeffs(i-2,j,k)& 02592 +ww1(2)*vv0+ww0(1)*vv1+ww1(1)*vv2 02593 coarse_coeffs(i-1,j,k)=coarse_coeffs(i-1,j,k)& 02594 +ww1(3)*vv0+ww0(2)*vv1+ww1(2)*vv2+ww0(1)*vv3+ww1(1)*vv4 02595 coarse_coeffs(i,j,k)=coarse_coeffs(i,j,k)& 02596 +ww1(4)*vv0+ww0(3)*vv1+ww1(3)*vv2+ww0(2)*vv3+ww1(2)*vv4+ww0(1)*vv5+ww1(1)*vv6 02597 coarse_coeffs(i+1,j,k)=coarse_coeffs(i+1,j,k)& 02598 +ww1(4)*vv2+ww0(3)*vv3+ww1(3)*vv4+ww0(2)*vv5+ww1(2)*vv6+ww0(1)*vv7 02599 coarse_coeffs(i+2,j,k)=coarse_coeffs(i+2,j,k)& 02600 +ww1(4)*vv4+ww0(3)*vv5+ww1(3)*vv6+ww0(2)*vv7 02601 coarse_coeffs(i+3,j,k)=coarse_coeffs(i+3,j,k)& 02602 +ww1(4)*vv6+ww0(3)*vv7 02603 END DO 02604 IF (.NOT.FLOOR((fine_bo(2,1)-f_shift(1))/2._dp)-coarse_bo(1,1)>=4) THEN 02605 CALL cp_assert(.FALSE.,& 02606 cp_failure_level,cp_assertion_failed,routineP,& 02607 "FLOOR((fine_bo(2,1)-f_shift(1))/2._dp)-coarse_bo(1,1)>=4",& 02608 error,failure) 02609 END IF 02610 rest_b=MODULO(FLOOR((fine_bo(2,1)-f_shift(1))/2._dp)-coarse_bo(1,1)-6,4) 02611 i=FLOOR((fine_bo(2,1)-f_shift(1))/2._dp)-3-rest_b+4 02612 CPPrecondition(fi==(i-2)*2+f_shift(1),cp_failure_level,routineP,error,failure) 02613 IF (rest_b>0) THEN 02614 fi=fi+1 02615 vv0=fine_values(fi,fj,fk) 02616 fi=fi+1 02617 vv1=fine_values(fi,fj,fk) 02618 coarse_coeffs(i-3,j,k)=coarse_coeffs(i-3,j,k)& 02619 +ww1(1)*vv0 02620 coarse_coeffs(i-2,j,k)=coarse_coeffs(i-2,j,k)& 02621 +ww1(2)*vv0+ww0(1)*vv1 02622 coarse_coeffs(i-1,j,k)=coarse_coeffs(i-1,j,k)& 02623 +ww1(3)*vv0+ww0(2)*vv1 02624 coarse_coeffs(i,j,k)=coarse_coeffs(i,j,k)& 02625 +ww1(4)*vv0+ww0(3)*vv1 02626 IF (rest_b>1) THEN 02627 fi=fi+1 02628 vv2=fine_values(fi,fj,fk) 02629 fi=fi+1 02630 vv3=fine_values(fi,fj,fk) 02631 coarse_coeffs(i-2,j,k)=coarse_coeffs(i-2,j,k)& 02632 +ww1(1)*vv2 02633 coarse_coeffs(i-1,j,k)=coarse_coeffs(i-1,j,k)& 02634 +ww1(2)*vv2+ww0(1)*vv3 02635 coarse_coeffs(i,j,k)=coarse_coeffs(i,j,k)& 02636 +ww1(3)*vv2+ww0(2)*vv3 02637 coarse_coeffs(i+1,j,k)=coarse_coeffs(i+1,j,k)& 02638 +ww1(4)*vv2+ww0(3)*vv3 02639 IF (rest_b>2) THEN 02640 fi=fi+1 02641 vv4=fine_values(fi,fj,fk) 02642 fi=fi+1 02643 vv5=fine_values(fi,fj,fk) 02644 fi=fi+1 02645 vv6=fine_values(fi,fj,fk) 02646 fi=fi+1 02647 vv7=fine_values(fi,fj,fk) 02648 coarse_coeffs(i-1,j,k)=coarse_coeffs(i-1,j,k)& 02649 +ww1(1)*vv4 02650 IF (has_i_ubound) THEN 02651 IF (coarse_bo(2,1)-2==FLOOR((fine_bo(2,1)-f_shift(1))/2._dp)) THEN 02652 fi=fi+1 02653 vv0=fine_values(fi,fj,fk) 02654 coarse_coeffs(i+4,j,k)=coarse_coeffs(i+4,j,k)& 02655 +vv0*ww1(4) 02656 ELSE 02657 vv0=0._dp 02658 END IF 02659 coarse_coeffs(i,j,k)=coarse_coeffs(i,j,k)& 02660 +ww1(2)*vv4+ww0(1)*vv5+ww1(1)*vv6 02661 coarse_coeffs(i+1,j,k)=coarse_coeffs(i+1,j,k)& 02662 +ww1(3)*vv4+ww0(2)*vv5+ww1(2)*vv6+ww0(1)*vv7+vv0*ww1(1) 02663 coarse_coeffs(i+2,j,k)=coarse_coeffs(i+2,j,k)& 02664 +ww1(4)*vv4+ww0(3)*vv5+ww1(3)*vv6+ww0(2)*vv7+vv0*ww1(2) 02665 coarse_coeffs(i+3,j,k)=coarse_coeffs(i+3,j,k)& 02666 +ww1(4)*vv6+ww0(3)*vv7+vv0*ww1(3) 02667 ELSEIF (pbc.and..not.is_split) THEN 02668 fi=fi+1 02669 vv0=fine_values(fi,fj,fk) 02670 vv1=fine_values(fine_bo(1,1),fj,fk) 02671 vv2=fine_values(fine_bo(1,1)+1,fj,fk) 02672 coarse_coeffs(i,j,k)=coarse_coeffs(i,j,k)& 02673 +ww1(2)*vv4+ww0(1)*vv5+ww1(1)*vv6 02674 coarse_coeffs(i+1,j,k)=coarse_coeffs(i+1,j,k)& 02675 +ww1(3)*vv4+ww0(2)*vv5+ww1(2)*vv6+ww0(1)*vv7+vv0*ww1(1) 02676 coarse_coeffs(i+2,j,k)=coarse_coeffs(i+2,j,k)& 02677 +ww1(4)*vv4+ww0(3)*vv5+ww1(3)*vv6+ww0(2)*vv7+vv0*ww1(2)& 02678 +vv1*ww0(1)+vv2*ww1(1) 02679 ELSE 02680 coarse_coeffs(i,j,k)=coarse_coeffs(i,j,k)& 02681 +ww1(2)*vv4+ww0(1)*vv5+wj*w_border1(3)*vv6 02682 coarse_coeffs(i+1,j,k)=coarse_coeffs(i+1,j,k)& 02683 +ww1(3)*vv4+ww0(2)*vv5+wj*w_border1(2)*vv6 02684 coarse_coeffs(i+2,j,k)=coarse_coeffs(i+2,j,k)& 02685 +ww1(4)*vv4+ww0(3)*vv5+wj*w_border1(1)*vv6+w_border0*wj*vv7 02686 END IF 02687 ELSE 02688 fi=fi+1 02689 vv4=fine_values(fi,fj,fk) 02690 fi=fi+1 02691 vv5=fine_values(fi,fj,fk) 02692 IF (has_i_ubound) THEN 02693 IF (coarse_bo(2,1)-2==FLOOR((fine_bo(2,1)-f_shift(1))/2._dp)) THEN 02694 fi=fi+1 02695 vv6=fine_values(fi,fj,fk) 02696 coarse_coeffs(i+3,j,k)=coarse_coeffs(i+3,j,k)& 02697 +ww1(4)*vv6 02698 ELSE 02699 vv6=0._dp 02700 END IF 02701 coarse_coeffs(i-1,j,k)=coarse_coeffs(i-1,j,k)& 02702 +ww1(1)*vv4 02703 coarse_coeffs(i,j,k)=coarse_coeffs(i,j,k)& 02704 +ww1(2)*vv4+ww0(1)*vv5+ww1(1)*vv6 02705 coarse_coeffs(i+1,j,k)=coarse_coeffs(i+1,j,k)& 02706 +ww1(3)*vv4+ww0(2)*vv5+ww1(2)*vv6 02707 coarse_coeffs(i+2,j,k)=coarse_coeffs(i+2,j,k)& 02708 +ww1(4)*vv4+ww0(3)*vv5+ww1(3)*vv6 02709 ELSEIF (pbc.AND..not.is_split) THEN 02710 fi=fi+1 02711 vv6=fine_values(fi,fj,fk) 02712 vv7=fine_values(fine_bo(1,1),fj,fk) 02713 vv0=fine_values(fine_bo(1,1)+1,fj,fk) 02714 coarse_coeffs(i-1,j,k)=coarse_coeffs(i-1,j,k)& 02715 +ww1(1)*vv4 02716 coarse_coeffs(i,j,k)=coarse_coeffs(i,j,k)& 02717 +ww1(4)*vv0+ww0(3)*vv1+ww1(3)*vv2+ww0(2)*vv3+ww1(2)*vv4+ww0(1)*vv5+ww1(1)*vv6 02718 coarse_coeffs(i+1,j,k)=coarse_coeffs(i+1,j,k)& 02719 +ww1(4)*vv2+ww0(3)*vv3+ww1(3)*vv4+ww0(2)*vv5+ww1(2)*vv6& 02720 +ww0(1)*vv7+ww1(1)*vv0 02721 ELSE 02722 coarse_coeffs(i-1,j,k)=coarse_coeffs(i-1,j,k)& 02723 +wj*w_border1(3)*vv4 02724 coarse_coeffs(i,j,k)=coarse_coeffs(i,j,k)& 02725 +wj*w_border1(2)*vv4 02726 coarse_coeffs(i+1,j,k)=coarse_coeffs(i+1,j,k)& 02727 +wj*(w_border1(1)*vv4+w_border0*vv5) 02728 END IF 02729 END IF 02730 ELSE 02731 fi=fi+1 02732 vv2=fine_values(fi,fj,fk) 02733 fi=fi+1 02734 vv3=fine_values(fi,fj,fk) 02735 IF (has_i_ubound) THEN 02736 IF (coarse_bo(2,1)-2==FLOOR((fine_bo(2,1)-f_shift(1))/2._dp)) THEN 02737 fi=fi+1 02738 vv4=fine_values(fi,fj,fk) 02739 coarse_coeffs(i+2,j,k)=coarse_coeffs(i+2,j,k)& 02740 +ww1(4)*vv4 02741 ELSE 02742 vv4=0._dp 02743 END IF 02744 coarse_coeffs(i-2,j,k)=coarse_coeffs(i-2,j,k)& 02745 +ww1(1)*vv2 02746 coarse_coeffs(i-1,j,k)=coarse_coeffs(i-1,j,k)& 02747 +ww1(2)*vv2+ww0(1)*vv3+ww1(1)*vv4 02748 coarse_coeffs(i,j,k)=coarse_coeffs(i,j,k)& 02749 +ww1(3)*vv2+ww0(2)*vv3+ww1(2)*vv4 02750 coarse_coeffs(i+1,j,k)=coarse_coeffs(i+1,j,k)& 02751 +ww1(4)*vv2+ww0(3)*vv3+ww1(3)*vv4 02752 ELSEIF (pbc.AND..not.is_split) THEN 02753 fi=fi+1 02754 vv4=fine_values(fi,fj,fk) 02755 vv5=fine_values(fine_bo(1,1),fj,fk) 02756 vv6=fine_values(fine_bo(1,1)+1,fj,fk) 02757 coarse_coeffs(i-2,j,k)=coarse_coeffs(i-2,j,k)& 02758 +ww1(1)*vv2 02759 coarse_coeffs(i-1,j,k)=coarse_coeffs(i-1,j,k)& 02760 +ww1(2)*vv2+ww0(1)*vv3+ww1(1)*vv4 02761 coarse_coeffs(i,j,k)=coarse_coeffs(i,j,k)& 02762 +ww1(3)*vv2+ww0(2)*vv3+ww1(2)*vv4+vv5*ww0(1)+ww1(1)*vv6 02763 ELSE 02764 coarse_coeffs(i-2,j,k)=coarse_coeffs(i-2,j,k)& 02765 +wj*w_border1(3)*vv2 02766 coarse_coeffs(i-1,j,k)=coarse_coeffs(i-1,j,k)& 02767 +wj*w_border1(2)*vv2 02768 coarse_coeffs(i,j,k)=coarse_coeffs(i,j,k)& 02769 +wj*(w_border1(1)*vv2+w_border0*vv3) 02770 END IF 02771 END IF 02772 ELSE 02773 fi=fi+1 02774 vv0=fine_values(fi,fj,fk) 02775 fi=fi+1 02776 vv1=fine_values(fi,fj,fk) 02777 IF (has_i_ubound) THEN 02778 IF (coarse_bo(2,1)-2==FLOOR((fine_bo(2,1)-f_shift(1))/2._dp)) THEN 02779 fi=fi+1 02780 vv2=fine_values(fi,fj,fk) 02781 coarse_coeffs(i+1,j,k)=coarse_coeffs(i+1,j,k)& 02782 +ww1(4)*vv2 02783 ELSE 02784 vv2=0._dp 02785 END IF 02786 coarse_coeffs(i-3,j,k)=coarse_coeffs(i-3,j,k)& 02787 +ww1(1)*vv0 02788 coarse_coeffs(i-2,j,k)=coarse_coeffs(i-2,j,k)& 02789 +ww1(2)*vv0+ww0(1)*vv1+ww1(1)*vv2 02790 coarse_coeffs(i-1,j,k)=coarse_coeffs(i-1,j,k)& 02791 +ww1(3)*vv0+ww0(2)*vv1+ww1(2)*vv2 02792 coarse_coeffs(i,j,k)=coarse_coeffs(i,j,k)& 02793 +ww1(4)*vv0+ww0(3)*vv1+ww1(3)*vv2 02794 ELSEIF (pbc.AND..not.is_split) THEN 02795 fi=fi+1 02796 vv2=fine_values(fi,fj,fk) 02797 vv3=fine_values(fine_bo(1,1),fk,fk) 02798 vv4=fine_values(fine_bo(1,1)+1,fk,fk) 02799 coarse_coeffs(i-3,j,k)=coarse_coeffs(i-3,j,k)& 02800 +ww1(1)*vv0 02801 coarse_coeffs(i-2,j,k)=coarse_coeffs(i-2,j,k)& 02802 +ww1(2)*vv0+ww0(1)*vv1+ww1(1)*vv2 02803 coarse_coeffs(i-1,j,k)=coarse_coeffs(i-1,j,k)& 02804 +ww1(3)*vv0+ww0(2)*vv1+ww1(2)*vv2+ww0(1)*vv3+ww1(1)*vv4 02805 ELSE 02806 coarse_coeffs(i-3,j,k)=coarse_coeffs(i-3,j,k)& 02807 +wj*w_border1(3)*vv0 02808 coarse_coeffs(i-2,j,k)=coarse_coeffs(i-2,j,k)& 02809 +wj*w_border1(2)*vv0 02810 coarse_coeffs(i-1,j,k)=coarse_coeffs(i-1,j,k)& 02811 +wj*(w_border1(1)*vv0+w_border0*vv1) 02812 END IF 02813 END IF 02814 CPPostcondition(fi==fine_bo(2,1),cp_failure_level,routineP,error,failure) 02815 END IF 02816 END DO 02817 END DO 02818 END DO 02819 END DO 02820 02821 ! *** parallel case 02822 IF (is_split) THEN 02823 CALL timeset(routineN//"_comm",handle2) 02824 coarse_slice_size=(coarse_bo(2,2)-coarse_bo(1,2)+1)*& 02825 (coarse_bo(2,3)-coarse_bo(1,3)+1) 02826 n_procs=coarse_coeffs_pw%pw_grid%para%group_size 02827 ALLOCATE(send_size(0:n_procs-1),send_offset(0:n_procs-1),& 02828 sent_size(0:n_procs-1),rcv_size(0:n_procs-1), & 02829 rcv_offset(0:n_procs-1), pp_lb(0:n_procs-1),& 02830 pp_ub(0:n_procs-1),real_rcv_size(0:n_procs-1),stat=stat) 02831 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 02832 02833 ! ** send size count 02834 02835 pos_of_x => coarse_coeffs_pw%pw_grid%para%pos_of_x 02836 send_size=0 02837 DO x=coarse_bo(1,1),coarse_bo(2,1) 02838 p=pos_of_x(coarse_gbo(1,1)+MODULO(x-coarse_gbo(1,1),s(1))) 02839 send_size(p)=send_size(p)+coarse_slice_size 02840 END DO 02841 02842 ! ** rcv size count 02843 02844 pos_of_x => fine_values_pw%pw_grid%para%pos_of_x 02845 p_old=pos_of_x(fine_gbo(1,1)) 02846 pp_lb=fine_gbo(2,1) 02847 pp_ub=fine_gbo(2,1)-1 02848 pp_lb(p_old)=fine_gbo(1,1) 02849 DO x=fine_gbo(1,1),fine_gbo(2,1) 02850 p=pos_of_x(x) 02851 IF (p/=p_old) THEN 02852 pp_ub(p_old)=x-1 02853 pp_lb(p)=x 02854 p_old=p 02855 END IF 02856 END DO 02857 pp_ub(p_old)=fine_gbo(2,1) 02858 02859 DO ip=0,n_procs-1 02860 IF (pp_lb(ip)<=pp_ub(ip)) THEN 02861 pp_lb(ip)=FLOOR(REAL(pp_lb(ip)-f_shift(1),dp)/2._dp)-1 02862 pp_ub(ip)=FLOOR(REAL(pp_ub(ip)+1-f_shift(1),dp)/2._dp)+1 02863 ELSE 02864 pp_lb(ip)=coarse_gbo(2,1) 02865 pp_ub(ip)=coarse_gbo(2,1)-1 02866 END IF 02867 IF (.NOT.is_split.OR..NOT.pbc) THEN 02868 pp_lb(ip)=MAX(pp_lb(ip),coarse_gbo(1,1)) 02869 pp_ub(ip)=MIN(pp_ub(ip),coarse_gbo(2,1)) 02870 END IF 02871 END DO 02872 02873 rcv_size=0 02874 DO ip=0,n_procs-1 02875 DO x=pp_lb(ip),coarse_gbo(1,1)-1 02876 x_att=coarse_gbo(1,1)+MODULO(x-coarse_gbo(1,1),s(1)) 02877 IF (x_att>=my_coarse_bo(1,1).AND.x_att<=my_coarse_bo(2,1)) THEN 02878 rcv_size(ip)=rcv_size(ip)+coarse_slice_size 02879 END IF 02880 END DO 02881 rcv_size(ip)=rcv_size(ip)+coarse_slice_size*& 02882 MAX(0,& 02883 MIN(pp_ub(ip),my_coarse_bo(2,1))-MAX(pp_lb(ip),my_coarse_bo(1,1))+1) 02884 DO x=coarse_gbo(2,1)+1,pp_ub(ip) 02885 x_att=coarse_gbo(1,1)+MODULO(x-coarse_gbo(1,1),s(1)) 02886 IF (x_att>=my_coarse_bo(1,1).AND.x_att<=my_coarse_bo(2,1)) THEN 02887 rcv_size(ip)=rcv_size(ip)+coarse_slice_size 02888 END IF 02889 END DO 02890 END DO 02891 02892 ! ** offsets & alloc send-rcv 02893 02894 send_tot_size=0 02895 DO ip=0,n_procs-1 02896 send_offset(ip)=send_tot_size 02897 send_tot_size=send_tot_size+send_size(ip) 02898 END DO 02899 CALL cp_assert(send_tot_size==& 02900 (coarse_bo(2,1)-coarse_bo(1,1)+1)*coarse_slice_size,& 02901 cp_failure_level,cp_assertion_failed,routineP,& 02902 "Error calculating send_tot_size "//& 02903 CPSourceFileRef,& 02904 error,failure) 02905 ALLOCATE(send_buf(0:send_tot_size-1),stat=stat) 02906 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 02907 02908 rcv_tot_size=0 02909 DO ip=0,n_procs-1 02910 rcv_offset(ip)=rcv_tot_size 02911 rcv_tot_size=rcv_tot_size+rcv_size(ip) 02912 END DO 02913 ALLOCATE(rcv_buf(0:rcv_tot_size-1),stat=stat) 02914 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 02915 02916 ! ** fill send buffer 02917 02918 pos_of_x => coarse_coeffs_pw%pw_grid%para%pos_of_x 02919 p_old=pos_of_x(coarse_gbo(1,1) & 02920 +MODULO(coarse_bo(1,1)-coarse_gbo(1,1),s(1))) 02921 sent_size=send_offset 02922 ss=coarse_bo(2,1)-coarse_bo(1,1)+1 02923 DO x=coarse_bo(1,1),coarse_bo(2,1) 02924 p=pos_of_x(coarse_gbo(1,1)+MODULO(x-coarse_gbo(1,1),s(1))) 02925 CALL dcopy(coarse_slice_size,& 02926 coarse_coeffs(x,coarse_bo(1,2),& 02927 coarse_bo(1,3)),ss,send_buf(sent_size(p)),1) 02928 sent_size(p)=sent_size(p)+coarse_slice_size 02929 END DO 02930 02931 CALL cp_assert(ALL(sent_size(0:n_procs-2)==send_offset(1:n_procs-1)),& 02932 cp_failure_level,cp_assertion_failed,routineP,& 02933 "error 1 filling send buffer "//& 02934 CPSourceFileRef,& 02935 error,failure) 02936 CALL cp_assert(sent_size(n_procs-1)==send_tot_size,& 02937 cp_failure_level,cp_assertion_failed,routineP,& 02938 "error 2 filling send buffer "//& 02939 CPSourceFileRef,& 02940 error,failure) 02941 02942 IF (local_data) THEN 02943 DEALLOCATE(coarse_coeffs,stat=stat) 02944 CPPostcondition(stat==0,cp_warning_level,routineP,error,failure) 02945 ELSE 02946 NULLIFY(coarse_coeffs) 02947 END IF 02948 02949 CPPostcondition(ALL(sent_size(:n_procs-2)==send_offset(1:)),cp_failure_level,routineP,error,failure) 02950 CPPostcondition(sent_size(n_procs-1)==send_tot_size,cp_failure_level,routineP,error,failure) 02951 ! test send/rcv sizes 02952 CALL mp_alltoall(send_size,real_rcv_size,1,coarse_coeffs_pw%pw_grid%para%group) 02953 02954 CPAssert(ALL(real_rcv_size==rcv_size),cp_failure_level,routineP,error,failure) 02955 ! all2all 02956 CALL mp_alltoall( sb=send_buf, scount=send_size, sdispl=send_offset,& 02957 rb=rcv_buf, rcount=rcv_size, rdispl=rcv_offset, & 02958 group=coarse_coeffs_pw%pw_grid%para%group ) 02959 02960 ! ** sum & reorder rcv buffer 02961 02962 sent_size=rcv_offset 02963 DO ip=0,n_procs-1 02964 02965 DO x=pp_lb(ip),coarse_gbo(1,1)-1 02966 x_att=coarse_gbo(1,1)+MODULO(x-coarse_gbo(1,1),s(1)) 02967 IF (x_att>=my_coarse_bo(1,1).AND.x_att<=my_coarse_bo(2,1)) THEN 02968 ii=sent_size(ip) 02969 DO k=coarse_bo(1,3),coarse_bo(2,3) 02970 DO j=coarse_bo(1,2),coarse_bo(2,2) 02971 coarse_coeffs_pw%cr3d(x_att,j,k)=coarse_coeffs_pw%cr3d(x_att,j,k)+rcv_buf(ii) 02972 ii=ii+1 02973 END DO 02974 END DO 02975 sent_size(ip)=ii 02976 END IF 02977 END DO 02978 02979 ii=sent_size(ip) 02980 DO x_att=MAX(pp_lb(ip),my_coarse_bo(1,1)),MIN(pp_ub(ip),my_coarse_bo(2,1)) 02981 DO k=coarse_bo(1,3),coarse_bo(2,3) 02982 DO j=coarse_bo(1,2),coarse_bo(2,2) 02983 coarse_coeffs_pw%cr3d(x_att,j,k)=coarse_coeffs_pw%cr3d(x_att,j,k)+rcv_buf(ii) 02984 ii=ii+1 02985 END DO 02986 END DO 02987 END DO 02988 sent_size(ip)=ii 02989 02990 DO x=coarse_gbo(2,1)+1,pp_ub(ip) 02991 x_att=coarse_gbo(1,1)+MODULO(x-coarse_gbo(1,1),s(1)) 02992 IF (x_att>=my_coarse_bo(1,1).AND.x_att<=my_coarse_bo(2,1)) THEN 02993 ii=sent_size(ip) 02994 DO k=coarse_bo(1,3),coarse_bo(2,3) 02995 DO j=coarse_bo(1,2),coarse_bo(2,2) 02996 coarse_coeffs_pw%cr3d(x_att,j,k)=coarse_coeffs_pw%cr3d(x_att,j,k)+rcv_buf(ii) 02997 ii=ii+1 02998 END DO 02999 END DO 03000 sent_size(ip)=ii 03001 END IF 03002 END DO 03003 03004 END DO 03005 03006 CALL cp_assert(ALL(sent_size(0:n_procs-2)==rcv_offset(1:n_procs-1)),& 03007 cp_failure_level,cp_assertion_failed,routineP,& 03008 "error 1 handling the rcv buffer "//& 03009 CPSourceFileRef,& 03010 error,failure) 03011 CALL cp_assert(sent_size(n_procs-1)==rcv_tot_size,& 03012 cp_failure_level,cp_assertion_failed,routineP,& 03013 "error 2 handling the rcv buffer "//& 03014 CPSourceFileRef,& 03015 error,failure) 03016 03017 ! dealloc 03018 DEALLOCATE(send_size,send_offset, rcv_size, rcv_offset, stat=stat) 03019 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 03020 DEALLOCATE(send_buf,rcv_buf,real_rcv_size, stat=stat) 03021 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 03022 DEALLOCATE(pp_ub,pp_lb,stat=stat) 03023 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 03024 CALL timestop(handle2) 03025 ELSE 03026 CPPostcondition(.NOT.local_data,cp_failure_level,routineP,error,failure) 03027 END IF 03028 END IF 03029 03030 CALL timestop(handle) 03031 END SUBROUTINE add_fine2coarse 03032 03033 ! ***************************************************************************** 03043 SUBROUTINE pw_spline_precond_create(preconditioner,precond_kind,& 03044 pool,pbc,transpose,error) 03045 TYPE(pw_spline_precond_type), POINTER :: preconditioner 03046 INTEGER, INTENT(in) :: precond_kind 03047 TYPE(pw_pool_type), POINTER :: pool 03048 LOGICAL, INTENT(in) :: pbc, transpose 03049 TYPE(cp_error_type), INTENT(inout) :: error 03050 03051 CHARACTER(len=*), PARAMETER :: routineN = 'pw_spline_precond_create', 03052 routineP = moduleN//':'//routineN 03053 03054 INTEGER :: stat 03055 LOGICAL :: failure 03056 03057 failure=.FALSE. 03058 03059 ALLOCATE(preconditioner,stat=stat) 03060 CPPostcondition(stat==0,cp_fatal_level,routineP,error,failure) 03061 IF (.NOT. failure) THEN 03062 last_precond_id=last_precond_id+1 03063 preconditioner%id_nr=last_precond_id 03064 preconditioner%ref_count=1 03065 preconditioner%kind=no_precond 03066 preconditioner%pool => pool 03067 preconditioner%pbc=pbc 03068 preconditioner%transpose=transpose 03069 CALL pw_pool_retain(pool,error=error) 03070 CALL pw_spline_precond_set_kind(preconditioner,precond_kind,error=error) 03071 END IF 03072 END SUBROUTINE pw_spline_precond_create 03073 03074 ! ***************************************************************************** 03082 SUBROUTINE pw_spline_precond_set_kind(preconditioner,precond_kind,pbc,& 03083 transpose,error) 03084 TYPE(pw_spline_precond_type), POINTER :: preconditioner 03085 INTEGER, INTENT(in) :: precond_kind 03086 LOGICAL, INTENT(in), OPTIONAL :: pbc, transpose 03087 TYPE(cp_error_type), INTENT(inout) :: error 03088 03089 CHARACTER(len=*), PARAMETER :: routineN = 'pw_spline_precond_set_kind', 03090 routineP = moduleN//':'//routineN 03091 03092 LOGICAL :: do_3d_coeff, failure 03093 REAL(kind=dp) :: s 03094 03095 failure=.FALSE. 03096 03097 CPPrecondition(ASSOCIATED(preconditioner),cp_failure_level,routineP,error,failure) 03098 CPPrecondition(preconditioner%ref_count>0,cp_failure_level,routineP,error,failure) 03099 IF (.NOT. failure) THEN 03100 IF (PRESENT(transpose)) preconditioner%transpose=transpose 03101 do_3d_coeff=.FALSE. 03102 preconditioner%kind=precond_kind 03103 IF (PRESENT(pbc)) preconditioner%pbc=pbc 03104 SELECT CASE(precond_kind) 03105 CASE(no_precond) 03106 CASE (precond_spl3_aint2) 03107 preconditioner%coeffs_1d=(/ -1.66_dp*0.25_dp,1.66_dp,-1.66_dp*0.25_dp /) 03108 preconditioner%sharpen=.FALSE. 03109 preconditioner%normalize=.FALSE. 03110 do_3d_coeff=.TRUE. 03111 CASE (precond_spl3_3) 03112 preconditioner%coeffs_1d(1)=-0.25_dp*1.6_dp 03113 preconditioner%coeffs_1d(2)=1.6_dp 03114 preconditioner%coeffs_1d(3)=-0.25_dp*1.6_dp 03115 preconditioner%sharpen=.FALSE. 03116 preconditioner%normalize=.FALSE. 03117 do_3d_coeff=.TRUE. 03118 CASE (precond_spl3_2) 03119 preconditioner%coeffs_1d(1)=-0.26_dp*1.76_dp 03120 preconditioner%coeffs_1d(2)=1.76_dp 03121 preconditioner%coeffs_1d(3)=-0.26_dp*1.76_dp 03122 preconditioner%sharpen=.FALSE. 03123 preconditioner%normalize=.FALSE. 03124 do_3d_coeff=.TRUE. 03125 CASE (precond_spl3_aint) 03126 preconditioner%coeffs_1d=spl3_1d_coeffs0 03127 preconditioner%sharpen=.TRUE. 03128 preconditioner%normalize=.TRUE. 03129 do_3d_coeff=.TRUE. 03130 CASE (precond_spl3_1) 03131 preconditioner%coeffs_1d(1)=0.5_dp/3._dp**(1._dp/3._dp) 03132 preconditioner%coeffs_1d(2)=4._dp/3._dp**(1._dp/3._dp) 03133 preconditioner%coeffs_1d(3)=0.5_dp/3._dp**(1._dp/3._dp) 03134 preconditioner%sharpen=.TRUE. 03135 preconditioner%normalize=.FALSE. 03136 do_3d_coeff=.TRUE. 03137 CASE default 03138 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 03139 END SELECT 03140 IF (do_3d_coeff) THEN 03141 s=1._dp 03142 IF (preconditioner%sharpen) s=-1._dp 03143 preconditioner%coeffs(1)=& 03144 s*preconditioner%coeffs_1d(2)*& 03145 preconditioner%coeffs_1d(2)*& 03146 preconditioner%coeffs_1d(2) 03147 preconditioner%coeffs(2)=& 03148 s*preconditioner%coeffs_1d(1)*& 03149 preconditioner%coeffs_1d(2)*& 03150 preconditioner%coeffs_1d(2) 03151 preconditioner%coeffs(3)=& 03152 s*preconditioner%coeffs_1d(1)*& 03153 preconditioner%coeffs_1d(1)*& 03154 preconditioner%coeffs_1d(2) 03155 preconditioner%coeffs(4)=& 03156 s*preconditioner%coeffs_1d(1)*& 03157 preconditioner%coeffs_1d(1)*& 03158 preconditioner%coeffs_1d(1) 03159 IF (preconditioner%sharpen) THEN 03160 IF (preconditioner%normalize) THEN 03161 preconditioner%coeffs(1)=2._dp+& 03162 preconditioner%coeffs(1) 03163 ELSE 03164 preconditioner%coeffs(1)=-preconditioner%coeffs(1) 03165 END IF 03166 END IF 03167 END IF 03168 END IF 03169 END SUBROUTINE pw_spline_precond_set_kind 03170 03171 ! ***************************************************************************** 03178 SUBROUTINE pw_spline_precond_retain(preconditioner,error) 03179 TYPE(pw_spline_precond_type), POINTER :: preconditioner 03180 TYPE(cp_error_type), INTENT(inout) :: error 03181 03182 CHARACTER(len=*), PARAMETER :: routineN = 'pw_spline_precond_retain', 03183 routineP = moduleN//':'//routineN 03184 03185 LOGICAL :: failure 03186 03187 failure=.FALSE. 03188 03189 CPPrecondition(ASSOCIATED(preconditioner),cp_failure_level,routineP,error,failure) 03190 IF (.NOT. failure) THEN 03191 CPPreconditionNoFail(preconditioner%ref_count>1,cp_failure_level,routineP,error) 03192 preconditioner%ref_count=preconditioner%ref_count+1 03193 END IF 03194 END SUBROUTINE pw_spline_precond_retain 03195 03196 ! ***************************************************************************** 03203 SUBROUTINE pw_spline_precond_release(preconditioner,error) 03204 TYPE(pw_spline_precond_type), POINTER :: preconditioner 03205 TYPE(cp_error_type), INTENT(inout) :: error 03206 03207 CHARACTER(len=*), PARAMETER :: routineN = 'pw_spline_precond_release', 03208 routineP = moduleN//':'//routineN 03209 03210 INTEGER :: stat 03211 LOGICAL :: failure 03212 03213 failure=.FALSE. 03214 03215 IF (ASSOCIATED(preconditioner)) THEN 03216 CPPreconditionNoFail(preconditioner%ref_count>0,cp_failure_level,routineP,error) 03217 preconditioner%ref_count=preconditioner%ref_count-1 03218 IF (preconditioner%ref_count==0) THEN 03219 CALL pw_pool_release(preconditioner%pool,error=error) 03220 DEALLOCATE(preconditioner,stat=stat) 03221 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 03222 END IF 03223 END IF 03224 END SUBROUTINE pw_spline_precond_release 03225 03226 ! ***************************************************************************** 03236 SUBROUTINE pw_spline_do_precond(preconditioner,in_v,out_v,error) 03237 TYPE(pw_spline_precond_type), POINTER :: preconditioner 03238 TYPE(pw_type), POINTER :: in_v, out_v 03239 TYPE(cp_error_type), INTENT(inout) :: error 03240 03241 CHARACTER(len=*), PARAMETER :: routineN = 'pw_spline_do_precond', 03242 routineP = moduleN//':'//routineN 03243 03244 LOGICAL :: failure 03245 03246 failure=.FALSE. 03247 03248 CPPrecondition(ASSOCIATED(preconditioner),cp_failure_level,routineP,error,failure) 03249 CPPrecondition(preconditioner%ref_count>0,cp_failure_level,routineP,error,failure) 03250 IF (.NOT. failure) THEN 03251 SELECT CASE(preconditioner%kind) 03252 CASE(no_precond) 03253 CALL pw_copy(in_v,out_v,error=error) 03254 CASE (precond_spl3_aint,precond_spl3_1) 03255 CALL pw_zero(out_v,error=error) 03256 IF (preconditioner%pbc) THEN 03257 CALL pw_nn_smear_r(pw_in=in_v,pw_out=out_v,& 03258 coeffs=preconditioner%coeffs,error=error) 03259 ELSE 03260 CALL pw_nn_compose_r_no_pbc(weights_1d=preconditioner%coeffs_1d,& 03261 pw_in=in_v, pw_out=out_v, sharpen=preconditioner%sharpen,& 03262 normalize=preconditioner%normalize,& 03263 transpose=preconditioner%transpose,error=error) 03264 END IF 03265 CASE(precond_spl3_3,precond_spl3_2,precond_spl3_aint2) 03266 CALL pw_zero(out_v,error=error) 03267 IF (preconditioner%pbc) THEN 03268 CALL pw_nn_smear_r(pw_in=in_v,pw_out=out_v,& 03269 coeffs=preconditioner%coeffs,error=error) 03270 ELSE 03271 CALL pw_nn_compose_r_no_pbc(weights_1d=preconditioner%coeffs_1d,& 03272 pw_in=in_v, pw_out=out_v, sharpen=preconditioner%sharpen,& 03273 normalize=preconditioner%normalize, smooth_boundary=.TRUE.,& 03274 transpose=preconditioner%transpose,error=error) 03275 END IF 03276 CASE default 03277 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 03278 END SELECT 03279 END IF 03280 END SUBROUTINE pw_spline_do_precond 03281 03282 ! ***************************************************************************** 03301 FUNCTION find_coeffs(values,coeffs,linOp,preconditioner, pool, & 03302 eps_r,eps_x,max_iter,error) RESULT(res) 03303 TYPE(pw_type), POINTER :: values, coeffs 03304 INTERFACE 03305 ! ***************************************************************************** 03306 SUBROUTINE linOp(pw_in,pw_out,error) 03307 USE pw_types, ONLY: pw_type 03308 USE cp_error_handling, ONLY: cp_error_type 03309 TYPE(pw_type), POINTER :: pw_in,pw_out 03310 TYPE(cp_error_type), INTENT(inout) :: error 03311 END SUBROUTINE linOp 03312 END INTERFACE 03313 TYPE(pw_spline_precond_type), POINTER :: preconditioner 03314 TYPE(pw_pool_type), POINTER :: pool 03315 REAL(kind=dp), INTENT(in) :: eps_r, eps_x 03316 INTEGER, INTENT(in) :: max_iter 03317 TYPE(cp_error_type), INTENT(inout) :: error 03318 LOGICAL :: res 03319 03320 CHARACTER(len=*), PARAMETER :: routineN = 'find_coeffs', 03321 routineP = moduleN//':'//routineN 03322 03323 INTEGER :: i, iiter, iter, j, k 03324 INTEGER, DIMENSION(2, 3) :: bo 03325 LOGICAL :: failure, last 03326 REAL(kind=dp) :: alpha, beta, eps_r_att, 03327 eps_x_att, r_z, r_z_new 03328 TYPE(cp_logger_type), POINTER :: logger 03329 TYPE(pw_type), POINTER :: Ap, p, r, z 03330 03331 failure=.FALSE. 03332 last=.FALSE. 03333 03334 res=.FALSE. 03335 IF (.NOT. failure) THEN 03336 logger => cp_error_get_logger(error) 03337 CALL pw_pool_create_pw(pool,r,use_data=REALDATA3D,in_space=REALSPACE,& 03338 error=error) 03339 CALL pw_pool_create_pw(pool,z,use_data=REALDATA3D,in_space=REALSPACE,& 03340 error=error) 03341 CALL pw_pool_create_pw(pool,p,use_data=REALDATA3D,in_space=REALSPACE,& 03342 error=error) 03343 CALL pw_pool_create_pw(pool,Ap,use_data=REALDATA3D,in_space=REALSPACE,& 03344 error=error) 03345 03346 CALL cp_add_iter_level(logger%iter_info,level_name="SPLINE_FIND_COEFFS",error=error) 03347 ext_do:DO iiter=1,max_iter,10 03348 CALL pw_zero(r,error=error) 03349 CALL linOp(pw_in=coeffs,pw_out=r,error=error) 03350 r%cr3d=-r%cr3d 03351 CALL pw_axpy(values,r,error=error) 03352 CALL pw_spline_do_precond(preconditioner,in_v=r,out_v=z,error=error) 03353 CALL pw_copy(z,p,error=error) 03354 r_z=pw_integral_ab_fast(r,z,error=error) 03355 03356 DO iter=iiter,MIN(iiter+9,max_iter) 03357 eps_r_att=SQRT(pw_integral_aa_fast(r,SQUARE,error=error)) 03358 IF (eps_r_att==0._dp) THEN 03359 eps_x_att=0._dp 03360 last=.TRUE. 03361 ELSE 03362 CALL pw_zero(Ap,error=error) 03363 CALL linOp(pw_in=p,pw_out=Ap,error=error) 03364 alpha=r_z/pw_integral_ab_fast(Ap,p,error=error) 03365 03366 CALL pw_axpy(p,coeffs,alpha=alpha,error=error) 03367 03368 eps_x_att=alpha*SQRT(pw_integral_aa_fast(p,SQUARE,error=error)) ! try to spare if unneded? 03369 IF (eps_r_att<eps_r .AND. eps_x_att<eps_x) last=.TRUE. 03370 END IF 03371 CALL cp_iterate(logger%iter_info,last=last,error=error) 03372 IF (last) THEN 03373 res=.TRUE. 03374 EXIT ext_do 03375 END IF 03376 03377 CALL pw_axpy(Ap,r,alpha=-alpha,error=error) 03378 03379 CALL pw_spline_do_precond(preconditioner,in_v=r,out_v=z,error=error) 03380 03381 r_z_new=pw_integral_ab_fast(r,z,error=error) 03382 beta=r_z_new/r_z 03383 r_z=r_z_new 03384 03385 bo=p%pw_grid%bounds_local 03386 DO k=bo(1,3),bo(2,3) 03387 DO j=bo(1,2),bo(2,2) 03388 DO i=bo(1,1),bo(2,1) 03389 p%cr3d(i,j,k)=z%cr3d(i,j,k)+beta*p%cr3d(i,j,k) 03390 END DO 03391 END DO 03392 END DO 03393 03394 END DO 03395 END DO ext_do 03396 CALL cp_rm_iter_level(logger%iter_info,level_name="SPLINE_FIND_COEFFS",error=error) 03397 03398 CALL pw_pool_give_back_pw(pool,r,error=error) 03399 CALL pw_pool_give_back_pw(pool,z,error=error) 03400 CALL pw_pool_give_back_pw(pool,p,error=error) 03401 CALL pw_pool_give_back_pw(pool,Ap,error=error) 03402 END IF 03403 END FUNCTION find_coeffs 03404 03405 ! ***************************************************************************** 03415 SUBROUTINE pw_restrict_s3(pw_fine_in,pw_coarse_out,coarse_pool,param_section,& 03416 error) 03417 TYPE(pw_type), POINTER :: pw_fine_in, pw_coarse_out 03418 TYPE(pw_pool_type), POINTER :: coarse_pool 03419 TYPE(section_vals_type), POINTER :: param_section 03420 TYPE(cp_error_type), INTENT(inout) :: error 03421 03422 CHARACTER(len=*), PARAMETER :: routineN = 'pw_restrict_s3', 03423 routineP = moduleN//':'//routineN 03424 03425 INTEGER :: aint_precond, handle, 03426 interp_kind, max_iter, 03427 precond_kind 03428 INTEGER, DIMENSION(2, 3) :: bo 03429 INTEGER, SAVE :: ifile = 0 03430 LOGICAL :: failure, pbc, 03431 safe_computation, success 03432 REAL(kind=dp) :: eps_r, eps_x 03433 TYPE(pw_spline_precond_type), POINTER :: precond 03434 TYPE(pw_type), POINTER :: coeffs, values 03435 03436 failure=.FALSE. 03437 ifile=ifile+1 03438 CALL timeset(routineN,handle) 03439 IF (.NOT. failure) THEN 03440 CALL section_vals_val_get(param_section,"safe_computation", & 03441 l_val=safe_computation, error=error) 03442 CALL section_vals_val_get(param_section,"aint_precond", & 03443 i_val=aint_precond, error=error) 03444 CALL section_vals_val_get(param_section,"precond", & 03445 i_val=precond_kind, error=error) 03446 CALL section_vals_val_get(param_section,"max_iter", & 03447 i_val=max_iter, error=error) 03448 CALL section_vals_val_get(param_section,"eps_r", & 03449 r_val=eps_r, error=error) 03450 CALL section_vals_val_get(param_section,"eps_x", & 03451 r_val=eps_x, error=error) 03452 CALL section_vals_val_get(param_section,"kind",& 03453 i_val=interp_kind, error=error) 03454 03455 pbc=(interp_kind==spline3_pbc_interp) 03456 CPPrecondition(pbc.OR.interp_kind==spline3_nopbc_interp,cp_failure_level,routineP,error,failure) 03457 bo=pw_coarse_out%pw_grid%bounds_local 03458 NULLIFY(values,coeffs) 03459 CALL pw_pool_create_pw(coarse_pool,values, use_data=REALDATA3D,& 03460 in_space=REALSPACE,error=error) 03461 CALL pw_zero(values,error=error) 03462 03463 !FM nullify(tst_pw) 03464 !FM CALL pw_pool_create_pw(coarse_pool,tst_pw, use_data=REALDATA3D,& 03465 !FM in_space=REALSPACE,error=error) 03466 !FM call pw_copy(values,tst_pw,error=error) 03467 !FM call add_fine2coarse(fine_values_pw=pw_fine_in,& 03468 !FM coarse_coeffs_pw=tst_pw,& 03469 !FM weights_1d=spl3_1d_transf_coeffs/2._dp, w_border0=0.5_dp,& 03470 !FM w_border1=spl3_1d_transf_border1/2._dp,pbc=pbc,& 03471 !FM safe_computation=.false.,error=error) 03472 03473 CALL add_fine2coarse(fine_values_pw=pw_fine_in,& 03474 coarse_coeffs_pw=values,& 03475 weights_1d=spl3_1d_transf_coeffs/2._dp, w_border0=0.5_dp,& 03476 w_border1=spl3_1d_transf_border1/2._dp,pbc=pbc,& 03477 safe_computation=safe_computation,error=error) 03478 03479 !FM CALL pw_compare_debug(tst_pw,values,max_diff,error=error) 03480 !FM WRITE(cp_logger_get_default_unit_nr(logger,.TRUE.),*)"f2cmax_diff=",max_diff 03481 !FM CALL pw_pool_give_back_pw(coarse_pool,tst_pw,error=error) 03482 03483 CALL pw_pool_create_pw(coarse_pool,coeffs, use_data=REALDATA3D,& 03484 in_space=REALSPACE,error=error) 03485 NULLIFY(precond) 03486 CALL pw_spline_precond_create(precond,precond_kind=aint_precond,& 03487 pool=coarse_pool,pbc=pbc,transpose=.TRUE.,error=error) 03488 CALL pw_spline_do_precond(precond,values,coeffs,error=error) 03489 CALL pw_spline_precond_set_kind(precond,precond_kind,error=error) 03490 IF (pbc) THEN 03491 success=find_coeffs(values=values,coeffs=coeffs,& 03492 linOp=spl3_pbc,preconditioner=precond, pool=coarse_pool, & 03493 eps_r=eps_r,eps_x=eps_x, max_iter=max_iter,error=error) 03494 ELSE 03495 success=find_coeffs(values=values,coeffs=coeffs,& 03496 linOp=spl3_nopbct,preconditioner=precond, pool=coarse_pool, & 03497 eps_r=eps_r,eps_x=eps_x, max_iter=max_iter,error=error) 03498 END IF 03499 CALL pw_spline_precond_release(precond,error=error) 03500 03501 CALL pw_zero(pw_coarse_out,error=error) 03502 CALL pw_axpy(coeffs,pw_coarse_out,error=error) 03503 03504 CALL pw_pool_give_back_pw(coarse_pool,values,error=error) 03505 CALL pw_pool_give_back_pw(coarse_pool,coeffs,error=error) 03506 END IF 03507 CALL timestop(handle) 03508 END SUBROUTINE pw_restrict_s3 03509 03510 ! ***************************************************************************** 03520 SUBROUTINE pw_prolongate_s3(pw_coarse_in,pw_fine_out,coarse_pool,& 03521 param_section,error) 03522 TYPE(pw_type), POINTER :: pw_coarse_in, pw_fine_out 03523 TYPE(pw_pool_type), POINTER :: coarse_pool 03524 TYPE(section_vals_type), POINTER :: param_section 03525 TYPE(cp_error_type), INTENT(inout) :: error 03526 03527 CHARACTER(len=*), PARAMETER :: routineN = 'pw_prolongate_s3', 03528 routineP = moduleN//':'//routineN 03529 03530 INTEGER :: aint_precond, handle, 03531 interp_kind, max_iter, 03532 precond_kind 03533 INTEGER, DIMENSION(2, 3) :: bo 03534 INTEGER, SAVE :: ifile = 0 03535 LOGICAL :: failure, pbc, 03536 safe_computation, success 03537 REAL(kind=dp) :: eps_r, eps_x 03538 TYPE(pw_spline_precond_type), POINTER :: precond 03539 TYPE(pw_type), POINTER :: coeffs 03540 03541 failure=.FALSE. 03542 03543 ifile=ifile+1 03544 CALL timeset(routineN,handle) 03545 IF (.NOT. failure) THEN 03546 NULLIFY(coeffs) 03547 CALL pw_pool_create_pw(coarse_pool,coeffs, use_data=REALDATA3D,& 03548 in_space=REALSPACE,error=error) 03549 bo=pw_coarse_in%pw_grid%bounds_local 03550 CALL section_vals_val_get(param_section,"safe_computation", & 03551 l_val=safe_computation, error=error) 03552 CALL section_vals_val_get(param_section,"aint_precond", & 03553 i_val=aint_precond, error=error) 03554 CALL section_vals_val_get(param_section,"precond", & 03555 i_val=precond_kind, error=error) 03556 CALL section_vals_val_get(param_section,"max_iter", & 03557 i_val=max_iter, error=error) 03558 CALL section_vals_val_get(param_section,"eps_r", & 03559 r_val=eps_r, error=error) 03560 CALL section_vals_val_get(param_section,"eps_x", & 03561 r_val=eps_x, error=error) 03562 CALL section_vals_val_get(param_section,"kind",& 03563 i_val=interp_kind,error=error) 03564 03565 pbc=(interp_kind==spline3_pbc_interp) 03566 CPPrecondition(pbc.OR.interp_kind==spline3_nopbc_interp,cp_failure_level,routineP,error,failure) 03567 NULLIFY(precond) 03568 CALL pw_spline_precond_create(precond,precond_kind=aint_precond,& 03569 pool=coarse_pool,pbc=pbc,transpose=.FALSE.,error=error) 03570 CALL pw_spline_do_precond(precond,pw_coarse_in,coeffs,error=error) 03571 CALL pw_spline_precond_set_kind(precond,precond_kind,error=error) 03572 IF (pbc) THEN 03573 success=find_coeffs(values=pw_coarse_in,coeffs=coeffs,& 03574 linOp=spl3_pbc,preconditioner=precond, pool=coarse_pool, & 03575 eps_r=eps_r,eps_x=eps_x,& 03576 max_iter=max_iter,error=error) 03577 ELSE 03578 success=find_coeffs(values=pw_coarse_in,coeffs=coeffs,& 03579 linOp=spl3_nopbc,preconditioner=precond, pool=coarse_pool, & 03580 eps_r=eps_r,eps_x=eps_x,& 03581 max_iter=max_iter,error=error) 03582 END IF 03583 CPPostconditionNoFail(success,cp_warning_level,routineP,error) 03584 CALL pw_spline_precond_release(precond,error=error) 03585 03586 !FM nullify(tst_pw) 03587 !FM call pw_create(tst_pw, pw_fine_out%pw_grid, use_data=REALDATA3D,& 03588 !FM in_space=REALSPACE, error=error) 03589 !FM call pw_copy(pw_fine_out,tst_pw,error=error) 03590 !FM CALL add_coarse2fine(coarse_coeffs_pw=coeffs,& 03591 !FM fine_values_pw=tst_pw,& 03592 !FM weights_1d=spl3_1d_transf_coeffs,& 03593 !FM w_border0=1._dp,& 03594 !FM w_border1=spl3_1d_transf_border1,& 03595 !FM pbc=pbc,safe_computation=.false.,& 03596 !FM error=error) 03597 03598 CALL add_coarse2fine(coarse_coeffs_pw=coeffs,& 03599 fine_values_pw=pw_fine_out,& 03600 weights_1d=spl3_1d_transf_coeffs,& 03601 w_border0=1._dp,& 03602 w_border1=spl3_1d_transf_border1,& 03603 pbc=pbc,safe_computation=safe_computation,& 03604 error=error) 03605 03606 !FM CALL pw_compare_debug(tst_pw,pw_fine_out,max_diff,error=error) 03607 !FM WRITE(cp_logger_get_default_unit_nr(logger,.TRUE.),*)"c2fmax_diff=",max_diff 03608 !FM CALL pw_release(tst_pw,error=error) 03609 03610 CALL pw_pool_give_back_pw(coarse_pool,coeffs,error=error) 03611 03612 END IF 03613 CALL timestop(handle) 03614 END SUBROUTINE pw_prolongate_s3 03615 03616 ! ***************************************************************************** 03625 SUBROUTINE pw_nn_compose_r_no_pbc(weights_1d,pw_in,pw_out,& 03626 sharpen,normalize,transpose,smooth_boundary,error) 03627 REAL(kind=dp), DIMENSION(-1:1) :: weights_1d 03628 TYPE(pw_type), POINTER :: pw_in, pw_out 03629 LOGICAL, INTENT(in), OPTIONAL :: sharpen, normalize, 03630 transpose, smooth_boundary 03631 TYPE(cp_error_type), INTENT(inout) :: error 03632 03633 CHARACTER(len=*), PARAMETER :: routineN = 'pw_nn_compose_r_no_pbc', 03634 routineP = moduleN//':'//routineN 03635 REAL(kind=dp), PARAMETER :: w_border = 0.0_dp 03636 03637 INTEGER :: first_index, i, j, jw, k, kw, 03638 last_index, myj, myk, n_els, 03639 stat 03640 INTEGER, DIMENSION(2, 3) :: bo, gbo 03641 INTEGER, DIMENSION(3) :: s 03642 LOGICAL :: failure, has_l_boundary, has_u_boundary, is_split, 03643 my_normalize, my_sharpen, my_smooth_boundary, my_transpose 03644 REAL(kind=dp) :: in_val_f, in_val_l, 03645 in_val_tmp, w_j, w_k 03646 REAL(kind=dp), DIMENSION(-1:1) :: w 03647 REAL(kind=dp), DIMENSION(:, :), POINTER :: l_boundary, tmp, u_boundary 03648 REAL(kind=dp), DIMENSION(:, :, :), 03649 POINTER :: in_val, out_val 03650 03651 bo=pw_in%pw_grid%bounds_local 03652 gbo=pw_in%pw_grid%bounds 03653 in_val => pw_in%cr3d 03654 out_val => pw_out%cr3d 03655 my_sharpen=.FALSE. 03656 IF (PRESENT(sharpen)) my_sharpen=sharpen 03657 my_normalize=.FALSE. 03658 IF (PRESENT(normalize)) my_normalize=normalize 03659 my_transpose=.FALSE. 03660 IF (PRESENT(transpose)) my_transpose=transpose 03661 my_smooth_boundary=.FALSE. 03662 IF (PRESENT(smooth_boundary)) my_smooth_boundary=smooth_boundary 03663 CPAssert(.not.my_normalize.OR.my_sharpen,cp_failure_level,routineP,error,failure) 03664 CPAssert(.NOT.my_smooth_boundary.OR..NOT.my_sharpen,cp_failure_level,routineP,error,failure) 03665 DO i=1,3 03666 s(i)=bo(2,i)-bo(1,i)+1 03667 END DO 03668 IF (ANY(s<1)) go to 666 03669 is_split= ANY(pw_in%pw_grid%bounds_local(:,1) /= & 03670 pw_in%pw_grid%bounds(:,1)) 03671 has_l_boundary=(gbo(1,1)==bo(1,1)) 03672 has_u_boundary= (gbo(2,1)==bo(2,1)) 03673 IF (is_split) THEN 03674 ALLOCATE(l_boundary(bo(1,2):bo(2,2),bo(1,3):bo(2,3)),& 03675 u_boundary(bo(1,2):bo(2,2),bo(1,3):bo(2,3)),& 03676 tmp(bo(1,2):bo(2,2),bo(1,3):bo(2,3)),stat=stat) 03677 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 03678 tmp(:,:)=pw_in%cr3d(bo(2,1),:,:) 03679 CALL mp_sendrecv(tmp,pw_in%pw_grid%para%pos_of_x(& 03680 gbo(1,1)+MODULO(bo(2,1)+1-gbo(1,1),gbo(2,1)-gbo(1,1)+1)),& 03681 l_boundary,pw_in%pw_grid%para%pos_of_x(& 03682 gbo(1,1)+MODULO(bo(1,1)-1-gbo(1,1),gbo(2,1)-gbo(1,1)+1)),& 03683 pw_in%pw_grid%para%group) 03684 tmp(:,:)=pw_in%cr3d(bo(1,1),:,:) 03685 CALL mp_sendrecv(tmp,pw_in%pw_grid%para%pos_of_x(& 03686 gbo(1,1)+MODULO(bo(1,1)-1-gbo(1,1),gbo(2,1)-gbo(1,1)+1)),& 03687 u_boundary,pw_in%pw_grid%para%pos_of_x(& 03688 gbo(1,1)+MODULO(bo(2,1)+1-gbo(1,1),gbo(2,1)-gbo(1,1)+1)),& 03689 pw_in%pw_grid%para%group) 03690 DEALLOCATE(tmp,stat=stat) 03691 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 03692 END IF 03693 03694 n_els=s(1) 03695 IF (has_l_boundary) THEN 03696 n_els=n_els-1 03697 first_index=bo(1,1)+1 03698 ELSE 03699 first_index=bo(1,1) 03700 END IF 03701 IF (has_u_boundary) THEN 03702 n_els=n_els-1 03703 last_index=bo(2,1)-1 03704 ELSE 03705 last_index=bo(2,1) 03706 END IF 03707 !#omp parallel do default(none) private(k,kw,myk,j,jw,myj,in_val_f,& 03708 !#omp in_val_l) shared(zderiv,yderiv,bo,in_val,out_val,s,l_boundary,& 03709 !#omp u_boundary,weights,is_split) 03710 DO k=bo(1,3),bo(2,3) 03711 DO kw=-1,1 03712 myk=k+kw 03713 IF (my_transpose) THEN 03714 IF (k>=gbo(2,3)-1.OR.k<=gbo(1,3)+1) THEN 03715 IF (k==gbo(2,3).OR.k==gbo(1,3)) THEN 03716 IF (myk<gbo(2,3).AND.myk>gbo(1,3)) THEN 03717 w_k=weights_1d(kw) 03718 IF (my_smooth_boundary) THEN 03719 w_k=weights_1d(kw)/weights_1d(0) 03720 END IF 03721 ELSE IF (kw==0) THEN 03722 w_k=1._dp 03723 ELSE 03724 CYCLE 03725 END IF 03726 ELSE 03727 IF (myk==gbo(2,3).OR.myk==gbo(1,3)) CYCLE 03728 w_k=weights_1d(kw) 03729 END IF 03730 ELSE 03731 w_k=weights_1d(kw) 03732 END IF 03733 ELSE 03734 IF (k>=gbo(2,3)-1.OR.k<=gbo(1,3)+1) THEN 03735 IF (k==gbo(2,3).OR.k==gbo(1,3)) THEN 03736 IF (kw/=0) CYCLE 03737 w_k=1._dp 03738 ELSE 03739 IF (my_smooth_boundary.AND.((k==gbo(1,3)+1.AND.myk==gbo(1,3)).OR.& 03740 (k==gbo(2,3)-1.AND.myk==gbo(2,3)))) THEN 03741 w_k=weights_1d(kw)/weights_1d(0) 03742 ELSE 03743 w_k=weights_1d(kw) 03744 END IF 03745 END IF 03746 ELSE 03747 w_k=weights_1d(kw) 03748 END IF 03749 END IF 03750 DO j=bo(1,2),bo(2,2) 03751 DO jw=-1,1 03752 myj=j+jw 03753 IF (j<gbo(2,2)-1.AND.j>gbo(1,2)+1) THEN 03754 w_j=w_k*weights_1d(jw) 03755 ELSE 03756 IF (my_transpose) THEN 03757 IF (j==gbo(2,2).OR.j==gbo(1,2)) THEN 03758 IF (myj<gbo(2,2).AND.myj>gbo(1,2)) THEN 03759 w_j=weights_1d(jw)*w_k 03760 IF (my_smooth_boundary) THEN 03761 w_j=weights_1d(jw)/weights_1d(0)*w_k 03762 END IF 03763 ELSE IF (jw==0) THEN 03764 w_j=w_k 03765 ELSE 03766 CYCLE 03767 END IF 03768 ELSE 03769 IF (myj==gbo(2,2).OR.myj==gbo(1,2)) CYCLE 03770 w_j=w_k*weights_1d(jw) 03771 END IF 03772 ELSE 03773 IF (j==gbo(2,2).OR.j==gbo(1,2)) THEN 03774 IF (jw/=0) CYCLE 03775 w_j=w_k 03776 ELSE IF (my_smooth_boundary.AND.((j==gbo(1,2)+1.AND.myj==gbo(1,2)).OR.& 03777 (j==gbo(2,2)-1.AND.myj==gbo(2,2)))) THEN 03778 w_j=w_k*weights_1d(jw)/weights_1d(0) 03779 ELSE 03780 w_j=w_k*weights_1d(jw) 03781 END IF 03782 END IF 03783 END IF 03784 03785 IF (has_l_boundary) THEN 03786 IF (my_transpose) THEN 03787 IF (s(1)==1) THEN 03788 CPAssert(.not.has_u_boundary,cp_failure_level,routineP,error,failure) 03789 in_val_tmp=u_boundary(myj,myk) 03790 ELSE 03791 in_val_tmp=in_val(bo(1,1)+1,myj,myk) 03792 END IF 03793 IF (my_sharpen) THEN 03794 IF (kw==0.AND.jw==0) THEN 03795 IF (my_normalize) THEN 03796 out_val(bo(1,1),j,k)=out_val(bo(1,1),j,k)+& 03797 (2.0_dp-w_j)*in_val(bo(1,1),myj,myk)-& 03798 in_val_tmp*weights_1d(1)*w_j 03799 ELSE 03800 out_val(bo(1,1),j,k)=out_val(bo(1,1),j,k)+& 03801 in_val(bo(1,1),myj,myk)*w_j-& 03802 in_val_tmp*weights_1d(1)*w_j 03803 END IF 03804 ELSE 03805 out_val(bo(1,1),j,k)=out_val(bo(1,1),j,k)-& 03806 in_val(bo(1,1),myj,myk)*w_j-& 03807 in_val_tmp*weights_1d(1)*w_j 03808 END IF 03809 ELSE IF (my_smooth_boundary) THEN 03810 out_val(bo(1,1),j,k)=out_val(bo(1,1),j,k)+& 03811 w_j*(in_val(bo(1,1),myj,myk)+& 03812 in_val_tmp*weights_1d(1)/weights_1d(0)) 03813 ELSE 03814 out_val(bo(1,1),j,k)=out_val(bo(1,1),j,k)+& 03815 w_j*(in_val(bo(1,1),myj,myk)+& 03816 in_val_tmp*weights_1d(1)) 03817 END IF 03818 in_val_f=0.0_dp 03819 ELSE 03820 in_val_f=in_val(bo(1,1),myj,myk) 03821 IF (my_sharpen) THEN 03822 IF (kw==0.AND.jw==0) THEN 03823 IF (my_normalize) THEN 03824 out_val(bo(1,1),j,k)=out_val(bo(1,1),j,k)+& 03825 (2.0_dp-w_j)*in_val_f 03826 ELSE 03827 out_val(bo(1,1),j,k)=out_val(bo(1,1),j,k)+& 03828 in_val_f*w_j 03829 END IF 03830 ELSE 03831 out_val(bo(1,1),j,k)=out_val(bo(1,1),j,k)-& 03832 in_val_f*w_j 03833 END IF 03834 ELSE 03835 out_val(bo(1,1),j,k)=out_val(bo(1,1),j,k)+& 03836 in_val_f*w_j 03837 END IF 03838 END IF 03839 ELSE 03840 in_val_f=l_boundary(myj,myk) 03841 END IF 03842 IF (has_u_boundary) THEN 03843 IF (my_transpose) THEN 03844 in_val_l=in_val(bo(2,1),myj,myk) 03845 IF (s(1)==1) THEN 03846 CPAssert(.not.has_l_boundary,cp_failure_level,routineP,error,failure) 03847 in_val_tmp=l_boundary(myj,myk) 03848 ELSE 03849 in_val_tmp=in_val(bo(2,1)-1,myj,myk) 03850 END IF 03851 IF (my_sharpen) THEN 03852 IF (kw==0.AND.jw==0) THEN 03853 IF (my_normalize) THEN 03854 out_val(bo(2,1),j,k)=out_val(bo(2,1),j,k)+& 03855 in_val_l*(2._dp-w_j)-& 03856 in_val_tmp*weights_1d(1)*w_j 03857 ELSE 03858 out_val(bo(2,1),j,k)=out_val(bo(2,1),j,k)+& 03859 in_val_l*w_j-& 03860 in_val_tmp*weights_1d(1)*w_j 03861 END IF 03862 ELSE 03863 out_val(bo(2,1),j,k)=out_val(bo(2,1),j,k)-& 03864 w_j*in_val_l-& 03865 in_val_tmp*weights_1d(1)*w_j 03866 END IF 03867 ELSE IF (my_smooth_boundary) THEN 03868 out_val(bo(2,1),j,k)=out_val(bo(2,1),j,k)+& 03869 w_j*(in_val_l+in_val_tmp*weights_1d(1)/weights_1d(0)) 03870 ELSE 03871 out_val(bo(2,1),j,k)=out_val(bo(2,1),j,k)+& 03872 w_j*(in_val_l+in_val_tmp*weights_1d(1)) 03873 END IF 03874 in_val_l=0._dp 03875 ELSE 03876 in_val_l=in_val(bo(2,1),myj,myk) 03877 IF (my_sharpen) THEN 03878 IF (kw==0.AND.jw==0) THEN 03879 IF (my_normalize) THEN 03880 out_val(bo(2,1),j,k)=out_val(bo(2,1),j,k)+& 03881 in_val_l*(2._dp-w_j) 03882 ELSE 03883 out_val(bo(2,1),j,k)=out_val(bo(2,1),j,k)+& 03884 in_val_l*w_j 03885 END IF 03886 ELSE 03887 out_val(bo(2,1),j,k)=out_val(bo(2,1),j,k)-& 03888 w_j*in_val_l 03889 END IF 03890 ELSE 03891 out_val(bo(2,1),j,k)=out_val(bo(2,1),j,k)+& 03892 w_j*in_val_l 03893 END IF 03894 END IF 03895 ELSE 03896 in_val_l=u_boundary(myj,myk) 03897 END IF 03898 IF (last_index>=first_index) THEN 03899 IF (my_transpose) THEN 03900 IF (bo(1,1)-1==gbo(1,1)) THEN 03901 in_val_f=0._dp 03902 ELSE IF (bo(2,1)+1==gbo(2,1)) THEN 03903 in_val_l=0._dp 03904 END IF 03905 END IF 03906 IF (my_sharpen) THEN 03907 w=-weights_1d*w_j 03908 IF (kw==0.AND.jw==0) THEN 03909 IF (my_normalize) THEN 03910 w(0)=w(0)+2._dp 03911 ELSE 03912 w(0)=-w(0) 03913 END IF 03914 END IF 03915 ELSE 03916 w=weights_1d*w_j 03917 END IF 03918 IF (my_smooth_boundary.AND.(.NOT.my_transpose)) THEN 03919 IF (gbo(1,1)+1>=bo(1,1).AND.& 03920 gbo(1,1)+1<=bo(2,1).AND.gbo(2,1)-gbo(1,1)>2) THEN 03921 IF (gbo(1,1)>=bo(1,1)) THEN 03922 out_val(gbo(1,1)+1,j,k)=out_val(gbo(1,1)+1,j,k)+& 03923 in_val(gbo(1,1),myj,myk)*w_j*weights_1d(-1)*& 03924 (1._dp/weights_1d(0)-1._dp) 03925 ELSE 03926 out_val(gbo(1,1)+1,j,k)=out_val(gbo(1,1)+1,j,k)+& 03927 l_boundary(myj,myk)*w_j*weights_1d(-1)*& 03928 (1._dp/weights_1d(0)-1._dp) 03929 END IF 03930 END IF 03931 END IF 03932 CALL pw_compose_stripe(weights=w,& 03933 in_val=in_val(first_index:last_index,myj,myk),& 03934 in_val_first=in_val_f,in_val_last=in_val_l,& 03935 out_val=out_val(first_index:last_index,j,k),& 03936 n_el=n_els) 03937 !FM call pw_compose_stripe2(weights=w,& 03938 !FM in_val=in_val,& 03939 !FM in_val_first=in_val_f,in_val_last=in_val_l,& 03940 !FM out_val=out_val,& 03941 !FM first_val=first_index,last_val=last_index,& 03942 !FM myj=myj,myk=myk,j=j,k=k) 03943 IF (my_smooth_boundary.AND.(.NOT.my_transpose))THEN 03944 IF (gbo(2,1)-1>=bo(1,1).AND.& 03945 gbo(2,1)-1<=bo(2,1).AND.gbo(2,1)-gbo(1,1)>2) THEN 03946 IF (gbo(2,1)<=bo(2,1)) THEN 03947 out_val(gbo(2,1)-1,j,k)=out_val(gbo(2,1)-1,j,k)+& 03948 in_val(gbo(2,1),myj,myk)*w_j*weights_1d(1)*& 03949 (1._dp/weights_1d(0)-1._dp) 03950 ELSE 03951 out_val(gbo(2,1)-1,j,k)=out_val(gbo(2,1)-1,j,k)+& 03952 u_boundary(myj,myk)*w_j*weights_1d(1)*& 03953 (1._dp/weights_1d(0)-1._dp) 03954 END IF 03955 END IF 03956 END IF 03957 03958 END IF 03959 END DO 03960 END DO 03961 END DO 03962 END DO 03963 03964 IF (is_split) THEN 03965 DEALLOCATE(l_boundary,u_boundary,stat=stat) 03966 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 03967 END IF 03968 666 CONTINUE 03969 END SUBROUTINE pw_nn_compose_r_no_pbc 03970 03971 ! ***************************************************************************** 03972 SUBROUTINE spl3_nopbc(pw_in,pw_out,error) 03973 TYPE(pw_type), POINTER :: pw_in, pw_out 03974 TYPE(cp_error_type), INTENT(inout) :: error 03975 03976 CALL pw_zero(pw_out,error=error) 03977 CALL pw_nn_compose_r_no_pbc(weights_1d=spl3_1d_coeffs0,pw_in=pw_in,& 03978 pw_out=pw_out,sharpen=.FALSE.,normalize=.FALSE.,error=error) 03979 END SUBROUTINE spl3_nopbc 03980 03981 ! ***************************************************************************** 03982 SUBROUTINE spl3_nopbct(pw_in,pw_out,error) 03983 TYPE(pw_type), POINTER :: pw_in, pw_out 03984 TYPE(cp_error_type), INTENT(inout) :: error 03985 03986 CALL pw_zero(pw_out,error=error) 03987 CALL pw_nn_compose_r_no_pbc(weights_1d=spl3_1d_coeffs0,pw_in=pw_in,& 03988 pw_out=pw_out,sharpen=.FALSE.,normalize=.FALSE.,transpose=.TRUE.,& 03989 error=error) 03990 END SUBROUTINE spl3_nopbct 03991 03992 ! ***************************************************************************** 03993 SUBROUTINE spl3_pbc(pw_in,pw_out,error) 03994 TYPE(pw_type), POINTER :: pw_in, pw_out 03995 TYPE(cp_error_type), INTENT(inout) :: error 03996 03997 CALL pw_zero(pw_out,error=error) 03998 CALL pw_nn_smear_r(pw_in,pw_out,coeffs=spline3_coeffs,& 03999 error=error) 04000 END SUBROUTINE spl3_pbc 04001 04002 ! ***************************************************************************** 04013 FUNCTION Eval_Interp_Spl3_pbc(vec,pw,error) RESULT(val) 04014 REAL(KIND=dp), DIMENSION(3), INTENT(in) :: vec 04015 TYPE(pw_type), POINTER :: pw 04016 TYPE(cp_error_type), INTENT(inout) :: error 04017 REAL(KIND=dp) :: val 04018 04019 CHARACTER(len=*), PARAMETER :: routineN = 'Eval_Interp_Spl3_pbc', 04020 routineP = moduleN//':'//routineN 04021 04022 INTEGER :: i, ivec(3), j, k, npts(3) 04023 INTEGER, DIMENSION(2, 3) :: bo, bo_l 04024 INTEGER, DIMENSION(4) :: ii, ij, ik 04025 LOGICAL :: failure, my_mpsum 04026 REAL(KIND=dp) :: a1, a2, a3, b1, b2, b3, c1, c2, c3, d1, d2, d3, dr1, 04027 dr2, dr3, e1, e2, e3, f1, f2, f3, g1, g2, g3, h1, h2, h3, p1, p2, p3, 04028 q1, q2, q3, r1, r2, r3, s1, s2, s3, s4, t1, t2, t3, t4, u1, u2, u3, v1, 04029 v2, v3, v4, xd1, xd2, xd3 04030 REAL(KIND=dp), DIMENSION(4, 4, 4) :: box 04031 REAL(KIND=dp), DIMENSION(:, :, :), 04032 POINTER :: grid 04033 04034 NULLIFY(grid) 04035 failure=.FALSE. 04036 IF (.NOT.failure) THEN 04037 my_mpsum = (pw%pw_grid%para%mode/=PW_MODE_LOCAL) 04038 npts = pw%pw_grid%npts 04039 ivec = FLOOR(vec/pw%pw_grid%dr) 04040 dr1 = pw%pw_grid%dr(1) 04041 dr2 = pw%pw_grid%dr(2) 04042 dr3 = pw%pw_grid%dr(3) 04043 04044 xd1 = (vec(1)/dr1)-REAL(ivec(1),kind=dp) 04045 xd2 = (vec(2)/dr2)-REAL(ivec(2),kind=dp) 04046 xd3 = (vec(3)/dr3)-REAL(ivec(3),kind=dp) 04047 grid => pw%cr3d(:,:,:) 04048 bo = pw%pw_grid%bounds 04049 bo_l = pw%pw_grid%bounds_local 04050 04051 ik(1) = MODULO(ivec(3)-1,npts(3)) + bo(1,3) 04052 ik(2) = MODULO(ivec(3) ,npts(3)) + bo(1,3) 04053 ik(3) = MODULO(ivec(3)+1,npts(3)) + bo(1,3) 04054 ik(4) = MODULO(ivec(3)+2,npts(3)) + bo(1,3) 04055 04056 ij(1) = MODULO(ivec(2)-1,npts(2)) + bo(1,2) 04057 ij(2) = MODULO(ivec(2) ,npts(2)) + bo(1,2) 04058 ij(3) = MODULO(ivec(2)+1,npts(2)) + bo(1,2) 04059 ij(4) = MODULO(ivec(2)+2,npts(2)) + bo(1,2) 04060 04061 ii(1) = MODULO(ivec(1)-1,npts(1)) + bo(1,1) 04062 ii(2) = MODULO(ivec(1) ,npts(1)) + bo(1,1) 04063 ii(3) = MODULO(ivec(1)+1,npts(1)) + bo(1,1) 04064 ii(4) = MODULO(ivec(1)+2,npts(1)) + bo(1,1) 04065 04066 DO k = 1,4 04067 DO j = 1,4 04068 DO i = 1,4 04069 IF(& 04070 ii(i) >= bo_l(1,1) .AND. & 04071 ii(i) <= bo_l(2,1) .AND. & 04072 ij(j) >= bo_l(1,2) .AND. & 04073 ij(j) <= bo_l(2,2) .AND. & 04074 ik(k) >= bo_l(1,3) .AND. & 04075 ik(k) <= bo_l(2,3) & 04076 ) THEN 04077 box(i,j,k) = grid(ii(i) + 1 - bo_l(1,1),& 04078 ij(j) + 1 - bo_l(1,2),& 04079 ik(k) + 1 - bo_l(1,3)) 04080 ELSE 04081 box(i,j,k) = 0.0_dp 04082 END IF 04083 END DO 04084 END DO 04085 END DO 04086 04087 a1 = 3.0_dp + xd1 04088 a2 = a1*a1 04089 a3 = a2*a1 04090 b1 = 2.0_dp + xd1 04091 b2 = b1*b1 04092 b3 = b2*b1 04093 c1 = 1.0_dp + xd1 04094 c2 = c1*c1 04095 c3 = c2*c1 04096 d1 = xd1 04097 d2 = d1*d1 04098 d3 = d2*d1 04099 e1 = 3.0_dp + xd2 04100 e2 = e1*e1 04101 e3 = e2*e1 04102 f1 = 2.0_dp + xd2 04103 f2 = f1*f1 04104 f3 = f2*f1 04105 g1 = 1.0_dp + xd2 04106 g2 = g1*g1 04107 g3 = g2*g1 04108 h1 = xd2 04109 h2 = h1*h1 04110 h3 = h2*h1 04111 p1 = 3.0_dp + xd3 04112 p2 = p1*p1 04113 p3 = p2*p1 04114 q1 = 2.0_dp + xd3 04115 q2 = q1*q1 04116 q3 = q2*q1 04117 r1 = 1.0_dp + xd3 04118 r2 = r1*r1 04119 r3 = r2*r1 04120 u1 = xd3 04121 u2 = u1*u1 04122 u3 = u2*u1 04123 04124 t1 = 1.0_dp/6.0_dp * (64.0_dp - 48.0_dp*a1 + 12.0_dp*a2 - a3) 04125 t2 = -22.0_dp/3.0_dp + 10.0_dp*b1 - 4.0_dp*b2 + 0.5_dp*b3 04126 t3 = 2.0_dp/3.0_dp - 2.0_dp*c1 + 2.0_dp*c2 - 0.5_dp*c3 04127 t4 = 1.0_dp/6.0_dp*d3 04128 s1 = 1.0_dp/6.0_dp * (64.0_dp - 48.0_dp*e1 + 12.0_dp*e2 - e3) 04129 s2 = -22.0_dp/3.0_dp + 10.0_dp*f1 - 4.0_dp*f2 + 0.5_dp*f3 04130 s3 = 2.0_dp/3.0_dp - 2.0_dp*g1 + 2.0_dp*g2 - 0.5_dp*g3 04131 s4 = 1.0_dp/6.0_dp*h3 04132 v1 = 1.0_dp/6.0_dp * (64.0_dp - 48.0_dp*p1 + 12.0_dp*p2 - p3) 04133 v2 = -22.0_dp/3.0_dp + 10.0_dp*q1 - 4.0_dp*q2 + 0.5_dp*q3 04134 v3 = 2.0_dp/3.0_dp - 2.0_dp*r1 + 2.0_dp*r2 - 0.5_dp*r3 04135 v4 = 1.0_dp/6.0_dp*u3 04136 04137 val = (( box(1,1,1) * t1 + box(2,1,1) * t2 + box(3,1,1) * t3 + box(4,1,1) * t4 ) * s1 +& 04138 ( box(1,2,1) * t1 + box(2,2,1) * t2 + box(3,2,1) * t3 + box(4,2,1) * t4 ) * s2 +& 04139 ( box(1,3,1) * t1 + box(2,3,1) * t2 + box(3,3,1) * t3 + box(4,3,1) * t4 ) * s3 +& 04140 ( box(1,4,1) * t1 + box(2,4,1) * t2 + box(3,4,1) * t3 + box(4,4,1) * t4 ) * s4 ) * v1 +& 04141 (( box(1,1,2) * t1 + box(2,1,2) * t2 + box(3,1,2) * t3 + box(4,1,2) * t4 ) * s1 +& 04142 ( box(1,2,2) * t1 + box(2,2,2) * t2 + box(3,2,2) * t3 + box(4,2,2) * t4 ) * s2 +& 04143 ( box(1,3,2) * t1 + box(2,3,2) * t2 + box(3,3,2) * t3 + box(4,3,2) * t4 ) * s3 +& 04144 ( box(1,4,2) * t1 + box(2,4,2) * t2 + box(3,4,2) * t3 + box(4,4,2) * t4 ) * s4 ) * v2 +& 04145 (( box(1,1,3) * t1 + box(2,1,3) * t2 + box(3,1,3) * t3 + box(4,1,3) * t4 ) * s1 +& 04146 ( box(1,2,3) * t1 + box(2,2,3) * t2 + box(3,2,3) * t3 + box(4,2,3) * t4 ) * s2 +& 04147 ( box(1,3,3) * t1 + box(2,3,3) * t2 + box(3,3,3) * t3 + box(4,3,3) * t4 ) * s3 +& 04148 ( box(1,4,3) * t1 + box(2,4,3) * t2 + box(3,4,3) * t3 + box(4,4,3) * t4 ) * s4 ) * v3 +& 04149 (( box(1,1,4) * t1 + box(2,1,4) * t2 + box(3,1,4) * t3 + box(4,1,4) * t4 ) * s1 +& 04150 ( box(1,2,4) * t1 + box(2,2,4) * t2 + box(3,2,4) * t3 + box(4,2,4) * t4 ) * s2 +& 04151 ( box(1,3,4) * t1 + box(2,3,4) * t2 + box(3,3,4) * t3 + box(4,3,4) * t4 ) * s3 +& 04152 ( box(1,4,4) * t1 + box(2,4,4) * t2 + box(3,4,4) * t3 + box(4,4,4) * t4 ) * s4 ) * v4 04153 04154 IF (my_mpsum) CALL mp_sum ( val, pw%pw_grid%para%group ) 04155 END IF 04156 04157 END FUNCTION Eval_Interp_Spl3_pbc 04158 04159 ! ***************************************************************************** 04170 FUNCTION Eval_d_Interp_Spl3_pbc(vec,pw,error) RESULT(val) 04171 REAL(KIND=dp), DIMENSION(3), INTENT(in) :: vec 04172 TYPE(pw_type), POINTER :: pw 04173 TYPE(cp_error_type), INTENT(inout) :: error 04174 REAL(KIND=dp) :: val(3) 04175 04176 CHARACTER(len=*), PARAMETER :: routineN = 'Eval_d_Interp_Spl3_pbc', 04177 routineP = moduleN//':'//routineN 04178 04179 INTEGER :: i, ivec(3), j, k, npts(3) 04180 INTEGER, DIMENSION(2, 3) :: bo, bo_l 04181 INTEGER, DIMENSION(4) :: ii, ij, ik 04182 LOGICAL :: failure, my_mpsum 04183 REAL(KIND=dp) :: a1, a2, a3, b1, b2, b3, c1, c2, c3, d1, d2, d3, dr1, 04184 dr1i, dr2, dr2i, dr3, dr3i, e1, e2, e3, f1, f2, f3, g1, g2, g3, h1, h2, 04185 h3, p1, p2, p3, q1, q2, q3, r1, r2, r3, s1, s1d, s1o, s2, s2d, s2o, s3, 04186 s3d, s3o, s4, s4d, s4o, t1, t1d, t1o, t2, t2d, t2o, t3, t3d, t3o, t4, 04187 t4d, t4o, u1, u2, u3, v1, v1d, v1o, v2, v2d, v2o, v3, v3d, v3o, v4, 04188 v4d, v4o, xd1, xd2, xd3 04189 REAL(KIND=dp), DIMENSION(4, 4, 4) :: box 04190 REAL(KIND=dp), DIMENSION(:, :, :), 04191 POINTER :: grid 04192 04193 NULLIFY(grid) 04194 failure=.FALSE. 04195 IF (.NOT.failure) THEN 04196 my_mpsum = (pw%pw_grid%para%mode/=PW_MODE_LOCAL) 04197 npts = pw%pw_grid%npts 04198 ivec = FLOOR(vec/pw%pw_grid%dr) 04199 dr1 = pw%pw_grid%dr(1) 04200 dr2 = pw%pw_grid%dr(2) 04201 dr3 = pw%pw_grid%dr(3) 04202 dr1i = 1.0_dp / dr1 04203 dr2i = 1.0_dp / dr2 04204 dr3i = 1.0_dp / dr3 04205 xd1 = (vec(1)/dr1)-REAL(ivec(1),kind=dp) 04206 xd2 = (vec(2)/dr2)-REAL(ivec(2),kind=dp) 04207 xd3 = (vec(3)/dr3)-REAL(ivec(3),kind=dp) 04208 grid => pw%cr3d(:,:,:) 04209 bo = pw%pw_grid%bounds 04210 bo_l = pw%pw_grid%bounds_local 04211 04212 ik(1) = MODULO(ivec(3)-1,npts(3)) + bo(1,3) 04213 ik(2) = MODULO(ivec(3) ,npts(3)) + bo(1,3) 04214 ik(3) = MODULO(ivec(3)+1,npts(3)) + bo(1,3) 04215 ik(4) = MODULO(ivec(3)+2,npts(3)) + bo(1,3) 04216 04217 ij(1) = MODULO(ivec(2)-1,npts(2)) + bo(1,2) 04218 ij(2) = MODULO(ivec(2) ,npts(2)) + bo(1,2) 04219 ij(3) = MODULO(ivec(2)+1,npts(2)) + bo(1,2) 04220 ij(4) = MODULO(ivec(2)+2,npts(2)) + bo(1,2) 04221 04222 ii(1) = MODULO(ivec(1)-1,npts(1)) + bo(1,1) 04223 ii(2) = MODULO(ivec(1) ,npts(1)) + bo(1,1) 04224 ii(3) = MODULO(ivec(1)+1,npts(1)) + bo(1,1) 04225 ii(4) = MODULO(ivec(1)+2,npts(1)) + bo(1,1) 04226 04227 DO k = 1,4 04228 DO j = 1,4 04229 DO i = 1,4 04230 IF(& 04231 ii(i) >= bo_l(1,1) .AND. & 04232 ii(i) <= bo_l(2,1) .AND. & 04233 ij(j) >= bo_l(1,2) .AND. & 04234 ij(j) <= bo_l(2,2) .AND. & 04235 ik(k) >= bo_l(1,3) .AND. & 04236 ik(k) <= bo_l(2,3) & 04237 ) THEN 04238 box(i,j,k) = grid(ii(i) + 1 - bo_l(1,1),& 04239 ij(j) + 1 - bo_l(1,2),& 04240 ik(k) + 1 - bo_l(1,3)) 04241 ELSE 04242 box(i,j,k) = 0.0_dp 04243 END IF 04244 END DO 04245 END DO 04246 END DO 04247 04248 a1 = 3.0_dp + xd1 04249 a2 = a1*a1 04250 a3 = a2*a1 04251 b1 = 2.0_dp + xd1 04252 b2 = b1*b1 04253 b3 = b2*b1 04254 c1 = 1.0_dp + xd1 04255 c2 = c1*c1 04256 c3 = c2*c1 04257 d1 = xd1 04258 d2 = d1*d1 04259 d3 = d2*d1 04260 e1 = 3.0_dp + xd2 04261 e2 = e1*e1 04262 e3 = e2*e1 04263 f1 = 2.0_dp + xd2 04264 f2 = f1*f1 04265 f3 = f2*f1 04266 g1 = 1.0_dp + xd2 04267 g2 = g1*g1 04268 g3 = g2*g1 04269 h1 = xd2 04270 h2 = h1*h1 04271 h3 = h2*h1 04272 p1 = 3.0_dp + xd3 04273 p2 = p1*p1 04274 p3 = p2*p1 04275 q1 = 2.0_dp + xd3 04276 q2 = q1*q1 04277 q3 = q2*q1 04278 r1 = 1.0_dp + xd3 04279 r2 = r1*r1 04280 r3 = r2*r1 04281 u1 = xd3 04282 u2 = u1*u1 04283 u3 = u2*u1 04284 04285 t1o = 1.0_dp/6.0_dp * (64.0_dp - 48.0_dp*a1 + 12.0_dp*a2 - a3) 04286 t2o = -22.0_dp/3.0_dp + 10.0_dp*b1 - 4.0_dp*b2 + 0.5_dp*b3 04287 t3o = 2.0_dp/3.0_dp - 2.0_dp*c1 + 2.0_dp*c2 - 0.5_dp*c3 04288 t4o = 1.0_dp/6.0_dp*d3 04289 s1o = 1.0_dp/6.0_dp * (64.0_dp - 48.0_dp*e1 + 12.0_dp*e2 - e3) 04290 s2o = -22.0_dp/3.0_dp + 10.0_dp*f1 - 4.0_dp*f2 + 0.5_dp*f3 04291 s3o = 2.0_dp/3.0_dp - 2.0_dp*g1 + 2.0_dp*g2 - 0.5_dp*g3 04292 s4o = 1.0_dp/6.0_dp*h3 04293 v1o = 1.0_dp/6.0_dp * (64.0_dp - 48.0_dp*p1 + 12.0_dp*p2 - p3) 04294 v2o = -22.0_dp/3.0_dp + 10.0_dp*q1 - 4.0_dp*q2 + 0.5_dp*q3 04295 v3o = 2.0_dp/3.0_dp - 2.0_dp*r1 + 2.0_dp*r2 - 0.5_dp*r3 04296 v4o = 1.0_dp/6.0_dp*u3 04297 04298 t1d = -8.0_dp + 4.0_dp*a1 - 0.5_dp*a2 04299 t2d = 10.0_dp - 8.0_dp*b1 + 1.5_dp*b2 04300 t3d = -2.0_dp + 4.0_dp*c1 - 1.5_dp*c2 04301 t4d = 0.5_dp*d2 04302 s1d = -8.0_dp + 4.0_dp*e1 - 0.5_dp*e2 04303 s2d = 10.0_dp - 8.0_dp*f1 + 1.5_dp*f2 04304 s3d = -2.0_dp + 4.0_dp*g1 - 1.5_dp*g2 04305 s4d = 0.5_dp*h2 04306 v1d = -8.0_dp + 4.0_dp*p1 - 0.5_dp*p2 04307 v2d = 10.0_dp - 8.0_dp*q1 + 1.5_dp*q2 04308 v3d = -2.0_dp + 4.0_dp*r1 - 1.5_dp*r2 04309 v4d = 0.5_dp*u2 04310 04311 t1 = t1d*dr1i 04312 t2 = t2d*dr1i 04313 t3 = t3d*dr1i 04314 t4 = t4d*dr1i 04315 s1 = s1o 04316 s2 = s2o 04317 s3 = s3o 04318 s4 = s4o 04319 v1 = v1o 04320 v2 = v2o 04321 v3 = v3o 04322 v4 = v4o 04323 val(1) = (( box(1,1,1) * t1 + box(2,1,1) * t2 + box(3,1,1) * t3 + box(4,1,1) * t4 ) * s1 +& 04324 ( box(1,2,1) * t1 + box(2,2,1) * t2 + box(3,2,1) * t3 + box(4,2,1) * t4 ) * s2 +& 04325 ( box(1,3,1) * t1 + box(2,3,1) * t2 + box(3,3,1) * t3 + box(4,3,1) * t4 ) * s3 +& 04326 ( box(1,4,1) * t1 + box(2,4,1) * t2 + box(3,4,1) * t3 + box(4,4,1) * t4 ) * s4 ) * v1 +& 04327 (( box(1,1,2) * t1 + box(2,1,2) * t2 + box(3,1,2) * t3 + box(4,1,2) * t4 ) * s1 +& 04328 ( box(1,2,2) * t1 + box(2,2,2) * t2 + box(3,2,2) * t3 + box(4,2,2) * t4 ) * s2 +& 04329 ( box(1,3,2) * t1 + box(2,3,2) * t2 + box(3,3,2) * t3 + box(4,3,2) * t4 ) * s3 +& 04330 ( box(1,4,2) * t1 + box(2,4,2) * t2 + box(3,4,2) * t3 + box(4,4,2) * t4 ) * s4 ) * v2 +& 04331 (( box(1,1,3) * t1 + box(2,1,3) * t2 + box(3,1,3) * t3 + box(4,1,3) * t4 ) * s1 +& 04332 ( box(1,2,3) * t1 + box(2,2,3) * t2 + box(3,2,3) * t3 + box(4,2,3) * t4 ) * s2 +& 04333 ( box(1,3,3) * t1 + box(2,3,3) * t2 + box(3,3,3) * t3 + box(4,3,3) * t4 ) * s3 +& 04334 ( box(1,4,3) * t1 + box(2,4,3) * t2 + box(3,4,3) * t3 + box(4,4,3) * t4 ) * s4 ) * v3 +& 04335 (( box(1,1,4) * t1 + box(2,1,4) * t2 + box(3,1,4) * t3 + box(4,1,4) * t4 ) * s1 +& 04336 ( box(1,2,4) * t1 + box(2,2,4) * t2 + box(3,2,4) * t3 + box(4,2,4) * t4 ) * s2 +& 04337 ( box(1,3,4) * t1 + box(2,3,4) * t2 + box(3,3,4) * t3 + box(4,3,4) * t4 ) * s3 +& 04338 ( box(1,4,4) * t1 + box(2,4,4) * t2 + box(3,4,4) * t3 + box(4,4,4) * t4 ) * s4 ) * v4 04339 04340 t1 = t1o 04341 t2 = t2o 04342 t3 = t3o 04343 t4 = t4o 04344 s1 = s1d*dr2i 04345 s2 = s2d*dr2i 04346 s3 = s3d*dr2i 04347 s4 = s4d*dr2i 04348 v1 = v1o 04349 v2 = v2o 04350 v3 = v3o 04351 v4 = v4o 04352 val(2) = (( box(1,1,1) * t1 + box(2,1,1) * t2 + box(3,1,1) * t3 + box(4,1,1) * t4 ) * s1 +& 04353 ( box(1,2,1) * t1 + box(2,2,1) * t2 + box(3,2,1) * t3 + box(4,2,1) * t4 ) * s2 +& 04354 ( box(1,3,1) * t1 + box(2,3,1) * t2 + box(3,3,1) * t3 + box(4,3,1) * t4 ) * s3 +& 04355 ( box(1,4,1) * t1 + box(2,4,1) * t2 + box(3,4,1) * t3 + box(4,4,1) * t4 ) * s4 ) * v1 +& 04356 (( box(1,1,2) * t1 + box(2,1,2) * t2 + box(3,1,2) * t3 + box(4,1,2) * t4 ) * s1 +& 04357 ( box(1,2,2) * t1 + box(2,2,2) * t2 + box(3,2,2) * t3 + box(4,2,2) * t4 ) * s2 +& 04358 ( box(1,3,2) * t1 + box(2,3,2) * t2 + box(3,3,2) * t3 + box(4,3,2) * t4 ) * s3 +& 04359 ( box(1,4,2) * t1 + box(2,4,2) * t2 + box(3,4,2) * t3 + box(4,4,2) * t4 ) * s4 ) * v2 +& 04360 (( box(1,1,3) * t1 + box(2,1,3) * t2 + box(3,1,3) * t3 + box(4,1,3) * t4 ) * s1 +& 04361 ( box(1,2,3) * t1 + box(2,2,3) * t2 + box(3,2,3) * t3 + box(4,2,3) * t4 ) * s2 +& 04362 ( box(1,3,3) * t1 + box(2,3,3) * t2 + box(3,3,3) * t3 + box(4,3,3) * t4 ) * s3 +& 04363 ( box(1,4,3) * t1 + box(2,4,3) * t2 + box(3,4,3) * t3 + box(4,4,3) * t4 ) * s4 ) * v3 +& 04364 (( box(1,1,4) * t1 + box(2,1,4) * t2 + box(3,1,4) * t3 + box(4,1,4) * t4 ) * s1 +& 04365 ( box(1,2,4) * t1 + box(2,2,4) * t2 + box(3,2,4) * t3 + box(4,2,4) * t4 ) * s2 +& 04366 ( box(1,3,4) * t1 + box(2,3,4) * t2 + box(3,3,4) * t3 + box(4,3,4) * t4 ) * s3 +& 04367 ( box(1,4,4) * t1 + box(2,4,4) * t2 + box(3,4,4) * t3 + box(4,4,4) * t4 ) * s4 ) * v4 04368 04369 t1 = t1o 04370 t2 = t2o 04371 t3 = t3o 04372 t4 = t4o 04373 s1 = s1o 04374 s2 = s2o 04375 s3 = s3o 04376 s4 = s4o 04377 v1 = v1d*dr3i 04378 v2 = v2d*dr3i 04379 v3 = v3d*dr3i 04380 v4 = v4d*dr3i 04381 val(3) = (( box(1,1,1) * t1 + box(2,1,1) * t2 + box(3,1,1) * t3 + box(4,1,1) * t4 ) * s1 +& 04382 ( box(1,2,1) * t1 + box(2,2,1) * t2 + box(3,2,1) * t3 + box(4,2,1) * t4 ) * s2 +& 04383 ( box(1,3,1) * t1 + box(2,3,1) * t2 + box(3,3,1) * t3 + box(4,3,1) * t4 ) * s3 +& 04384 ( box(1,4,1) * t1 + box(2,4,1) * t2 + box(3,4,1) * t3 + box(4,4,1) * t4 ) * s4 ) * v1 +& 04385 (( box(1,1,2) * t1 + box(2,1,2) * t2 + box(3,1,2) * t3 + box(4,1,2) * t4 ) * s1 +& 04386 ( box(1,2,2) * t1 + box(2,2,2) * t2 + box(3,2,2) * t3 + box(4,2,2) * t4 ) * s2 +& 04387 ( box(1,3,2) * t1 + box(2,3,2) * t2 + box(3,3,2) * t3 + box(4,3,2) * t4 ) * s3 +& 04388 ( box(1,4,2) * t1 + box(2,4,2) * t2 + box(3,4,2) * t3 + box(4,4,2) * t4 ) * s4 ) * v2 +& 04389 (( box(1,1,3) * t1 + box(2,1,3) * t2 + box(3,1,3) * t3 + box(4,1,3) * t4 ) * s1 +& 04390 ( box(1,2,3) * t1 + box(2,2,3) * t2 + box(3,2,3) * t3 + box(4,2,3) * t4 ) * s2 +& 04391 ( box(1,3,3) * t1 + box(2,3,3) * t2 + box(3,3,3) * t3 + box(4,3,3) * t4 ) * s3 +& 04392 ( box(1,4,3) * t1 + box(2,4,3) * t2 + box(3,4,3) * t3 + box(4,4,3) * t4 ) * s4 ) * v3 +& 04393 (( box(1,1,4) * t1 + box(2,1,4) * t2 + box(3,1,4) * t3 + box(4,1,4) * t4 ) * s1 +& 04394 ( box(1,2,4) * t1 + box(2,2,4) * t2 + box(3,2,4) * t3 + box(4,2,4) * t4 ) * s2 +& 04395 ( box(1,3,4) * t1 + box(2,3,4) * t2 + box(3,3,4) * t3 + box(4,3,4) * t4 ) * s3 +& 04396 ( box(1,4,4) * t1 + box(2,4,4) * t2 + box(3,4,4) * t3 + box(4,4,4) * t4 ) * s4 ) * v4 04397 04398 IF (my_mpsum) CALL mp_sum ( val, pw%pw_grid%para%group ) 04399 04400 END IF 04401 04402 END FUNCTION Eval_d_Interp_Spl3_pbc 04403 04404 END MODULE pw_spline_utils 04405
1.7.3