Direct M2M translation over OZ axis
Compute the following matrix-vector product: \f[ \mathrm{dst} = \beta \mathrm{dst} + \alpha M_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$ M_M \f$ is a matrix of multipole-to-multipole translation over OZ axis.
@param[in] z: OZ coordinate from new to old centers of harmonics
@param[in] src_r: Radius of old harmonics
@param[in] dst_r: Radius of new harmonics
@parma[in] p: Maximal degree of spherical harmonics
@param[in] vscales: Normalization constants for harmonics
@param[in] vcnk: Square roots of combinatorial numbers C_n^k
@param[in] alpha: Scalar multipler for alpha
@param[in] src_m: Expansion in old harmonics
@param[in] beta: Scalar multipler for dst_m
@param[inout] dst_m: Expansion in new harmonics
@param[out] work: Temporary workspace of a size (2*(p+1))
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=rp), | intent(in) | :: | z | |||
| integer, | intent(in) | :: | p | |||
| real(kind=rp), | intent(in) | :: | vscales((p+1)*(p+1)) | |||
| real(kind=rp), | intent(in) | :: | vcnk((2*p+1)*(p+1)) | |||
| real(kind=rp), | intent(in) | :: | alpha | |||
| real(kind=rp), | intent(in) | :: | src_m((p+1)*(p+1)) | |||
| real(kind=rp), | intent(in) | :: | beta | |||
| real(kind=rp), | intent(inout) | :: | dst_m((p+1)*(p+1)) | |||
| real(kind=rp), | intent(out), | target | :: | work(2*(p+1)) |
subroutine fmm_m2m_ztranslate_work(z, p, vscales, vcnk, alpha, &
& src_m, beta, dst_m, work)
implicit none
! Inputs
integer, intent(in) :: p
real(rp), intent(in) :: z, vscales((p+1)*(p+1)), &
& vcnk((2*p+1)*(p+1)), alpha, src_m((p+1)*(p+1)), beta
! Output
real(rp), intent(inout) :: dst_m((p+1)*(p+1))
! Temporary workspace
real(rp), intent(out), target :: work(2*(p+1))
! Local variables
real(rp) :: r2, tmp1, tmp2, tmp3, res1, res2
integer :: j, k, n, indj, indjn, indjk1, indjk2
! Pointers for temporary values of powers
real(rp), pointer :: pow_r2(:)
! In case alpha is 0.0 just do a proper scaling of output
if (abs(alpha) < eps_rp) then
if (abs(beta) < eps_rp) then
dst_m = 0.0
else
dst_m = beta * dst_m
end if
return
end if
! Now alpha is non-0.0
! If harmonics have different centers
if (abs(z) > eps_rp) then
! Prepare pointers
pow_r2(1:p+1) => work(1:p+1)
! Get ratios r1 and r2
r2 = z
pow_r2(1) = 1.0
do j = 2, p+1
pow_r2(j) = pow_r2(j-1) * r2
end do
! Do actual M2M
! Overwrite output if beta is 0.0
if (abs(beta) < eps_rp) then
do j = 0, p
! Offset for dst_m
indj = j*j + j + 1
! k = 0
tmp1 = alpha * vscales(indj)
tmp2 = tmp1
res1 = 0.0
! Offset for vcnk
indjk1 = j*(j+1)/2 + 1
do n = 0, j
! Offset for src_m
indjn = (j-n)**2 + (j-n) + 1
tmp3 = pow_r2(n+1) / &
& vscales(indjn) * vcnk(indjk1+n)**2
res1 = res1 + tmp3*src_m(indjn)
end do
dst_m(indj) = tmp2 * res1
! k != 0
do k = 1, j
tmp2 = tmp1
res1 = 0.0
res2 = 0.0
! Offsets for vcnk
indjk1 = (j-k)*(j-k+1)/2 + 1
indjk2 = (j+k)*(j+k+1)/2 + 1
do n = 0, j-k
! Offset for src_m
indjn = (j-n)**2 + (j-n) + 1
tmp3 = pow_r2(n+1) / &
& vscales(indjn) * vcnk(indjk1+n) * &
& vcnk(indjk2+n)
res1 = res1 + tmp3*src_m(indjn+k)
res2 = res2 + tmp3*src_m(indjn-k)
end do
dst_m(indj+k) = tmp2 * res1
dst_m(indj-k) = tmp2 * res2
end do
end do
! Update output if beta is non-0.0
else
do j = 0, p
! Offset for dst_m
indj = j*j + j + 1
! k = 0
tmp1 = alpha * vscales(indj)
tmp2 = tmp1
res1 = 0.0
! Offset for vcnk
indjk1 = j * (j+1) /2 + 1
do n = 0, j
! Offset for src_m
indjn = (j-n)**2 + (j-n) + 1
tmp3 = pow_r2(n+1) / &
& vscales(indjn) * vcnk(indjk1+n)**2
res1 = res1 + tmp3*src_m(indjn)
end do
dst_m(indj) = beta*dst_m(indj) + tmp2*res1
! k != 0
do k = 1, j
tmp2 = tmp1
res1 = 0.0
res2 = 0.0
! Offsets for vcnk
indjk1 = (j-k)*(j-k+1)/2 + 1
indjk2 = (j+k)*(j+k+1)/2 + 1
do n = 0, j-k
! Offset for src_m
indjn = (j-n)**2 + (j-n) + 1
tmp3 = pow_r2(n+1) / &
& vscales(indjn) * vcnk(indjk1+n) * &
& vcnk(indjk2+n)
res1 = res1 + tmp3*src_m(indjn+k)
res2 = res2 + tmp3*src_m(indjn-k)
end do
dst_m(indj+k) = beta*dst_m(indj+k) + tmp2*res1
dst_m(indj-k) = beta*dst_m(indj-k) + tmp2*res2
end do
end do
end if
! If harmonics are located at the same point
else
! Overwrite output if beta is 0.0
if (abs(beta) < eps_rp) then
tmp1 = alpha
do j = 0, p
indj = j*j + j + 1
do k = indj-j, indj+j
dst_m(k) = tmp1 * src_m(k)
end do
end do
! Update output if beta is non-0.0
else
tmp1 = alpha
do j = 0, p
indj = j*j + j + 1
do k = indj-j, indj+j
dst_m(k) = beta*dst_m(k) + tmp1*src_m(k)
end do
tmp1 = tmp1
end do
end if
end if
end subroutine fmm_m2m_ztranslate_work