Direct M2L translation by 4 rotations and 1 translation
Compute the following matrix-vector product: \f[ \mathrm{dst} = \beta \mathrm{dst} + \alpha L_M \mathrm{src}, \f] where \f$ \mathrm{dst} \f$ is a vector of coefficients of output spherical harmonics, \f$ \mathrm{src} \f$ is a vector of coefficients of input spherical harmonics and \f$ L_M \f$ is a matrix of a multipole-to-local translation.
Rotates around OZ and OY axes, translates over OZ and then rotates back around OY and OZ axes.
@param[in] c: Radius-vector from new to old centers of harmonics
@param[in] src_r: Radius of old harmonics
@param[in] dst_r: Radius of new harmonics
@param[in] pm: Maximal degree of multipole spherical harmonics
@param[in] pl: Maximal degree of local spherical harmonics
@param[in] vscales: Normalization constants for Y_lm
@param[in] m2l_ztranslate_coef:
@param[in] alpha: Scalar multiplier for src_m
@param[in] src_m: Expansion in old harmonics
@param[in] beta: Scalar multiplier for dst_l
@param[inout] dst_l: Expansion in new harmonics
@param[out] work: Temporary workspace of a size 6pp+19*p+8 where p is a
maximum of pm and pl
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=rp), | intent(in) | :: | c(3) | |||
integer, | intent(in) | :: | pm | |||
integer, | intent(in) | :: | pl | |||
real(kind=rp), | intent(in) | :: | vscales((pm+pl+1)**2) | |||
real(kind=rp), | intent(in) | :: | m2l_ztranslate_coef(pm+1,pl+1,pl+1) | |||
real(kind=rp), | intent(in) | :: | alpha | |||
real(kind=rp), | intent(in) | :: | src_m((pm+1)*(pm+1)) | |||
real(kind=rp), | intent(in) | :: | beta | |||
real(kind=rp), | intent(inout) | :: | dst_l((pl+1)*(pl+1)) | |||
real(kind=rp), | intent(out), | target | :: | work(6*max(pm,pl)**2+19*max(pm,pl)+8) |
subroutine fmm_m2l_rotation_work(c, pm, pl, vscales, &
& m2l_ztranslate_coef, alpha, src_m, beta, dst_l, work)
! Inputs
integer, intent(in) :: pm, pl
real(rp), intent(in) :: c(3), vscales((pm+pl+1)**2), &
& m2l_ztranslate_coef(pm+1, pl+1, pl+1), alpha, src_m((pm+1)*(pm+1)), &
& beta
! Output
real(rp), intent(inout) :: dst_l((pl+1)*(pl+1))
! Temporary workspace
real(rp), intent(out), target :: &
& work(6*max(pm, pl)**2 + 19*max(pm, pl) + 8)
! Local variables
real(rp) :: rho, ctheta, stheta, cphi, sphi
integer :: m, n, p
! Pointers for temporary values of harmonics
real(rp), pointer :: tmp_ml(:), tmp_ml2(:), vcos(:), vsin(:)
! Covert Cartesian coordinates into spherical
call carttosph(c, rho, ctheta, stheta, cphi, sphi)
! If no need for rotations, just do translation along z
if (abs(stheta) < eps_rp) then
! Workspace here is (pm+2)*(pm+1)
call fmm_m2l_ztranslate_work(c(3), pm, pl, vscales, &
& m2l_ztranslate_coef, alpha, src_m, beta, dst_l, work)
return
end if
! Prepare pointers
p = max(pm, pl)
m = (p+1)**2
n = 4*m + 5*p ! 4*p*p + 13*p + 4
tmp_ml(1:m) => work(n+1:n+m) ! 5*p*p + 15*p + 5
n = n + m
tmp_ml2(1:m) => work(n+1:n+m) ! 6*p*p + 17*p + 6
n = n + m
m = p + 1
vcos => work(n+1:n+m) ! 6*p*p + 18*p + 7
n = n + m
vsin => work(n+1:n+m) ! 6*p*p + 19*p + 8
! Compute arrays of cos and sin that are needed for rotations of harmonics
call trgev(cphi, sphi, p, vcos, vsin)
! Rotate around OZ axis (work array might appear in the future)
call fmm_sph_rotate_oz_adj_work(pm, vcos, vsin, alpha, src_m, 0.0_rp, tmp_ml)
! Perform rotation in the OXZ plane, work size is 4*pm*pm+13*pm+4
call fmm_sph_rotate_oxz_work(pm, ctheta, -stheta, 1.0_rp, tmp_ml, 0.0_rp, &
& tmp_ml2, work)
! OZ translation, workspace here is (pm+2)*(pm+1)
call fmm_m2l_ztranslate_work(rho, pm, pl, vscales, &
& m2l_ztranslate_coef, 1.0_rp, tmp_ml2, 0.0_rp, tmp_ml, work)
! Backward rotation in the OXZ plane, work size is 4*pl*pl+13*pl+4
call fmm_sph_rotate_oxz_work(pl, ctheta, stheta, 1.0_rp, tmp_ml, 0.0_rp, &
& tmp_ml2, work)
! Backward rotation around OZ axis (work array might appear in the future)
call fmm_sph_rotate_oz_work(pl, vcos, vsin, 1.0_rp, tmp_ml2, beta, dst_l)
end subroutine fmm_m2l_rotation_work