fmm_l2l_ztranslate_work Subroutine

private subroutine fmm_l2l_ztranslate_work(z, p, vscales, vfact, alpha, src_l, beta, dst_l, work)

Direct L2L translation over OZ axis

Compute the following matrix-vector product: \f[ \mathrm{dst} = \beta \mathrm{dst} + \alpha L_L \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_L \f$ is a matrix of local-to-local 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] vfact: Square roots of factorials @param[in] alpha: Scalar multipler for src_l @param[in] src_l: Expansion in old harmonics @param[in] beta: Scalar multipler for dst_l @param[inout] dst_l: Expansion in new harmonics @param[out] work: Temporary workspace of a size (2*(p+1))

Arguments

Type IntentOptional 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) :: vfact(2*p+1)
real(kind=rp), intent(in) :: alpha
real(kind=rp), intent(in) :: src_l((p+1)*(p+1))
real(kind=rp), intent(in) :: beta
real(kind=rp), intent(inout) :: dst_l((p+1)*(p+1))
real(kind=rp), intent(out), target :: work(2*(p+1))

Called by

proc~~fmm_l2l_ztranslate_work~~CalledByGraph proc~fmm_l2l_ztranslate_work fmm_l2l_ztranslate_work proc~fmm_l2l_rotation_work fmm_l2l_rotation_work proc~fmm_l2l_rotation_work->proc~fmm_l2l_ztranslate_work proc~fmm_l2l fmm_l2l proc~fmm_l2l->proc~fmm_l2l_rotation_work proc~tree_l2l tree_l2l proc~tree_l2l->proc~fmm_l2l proc~cart_propfar_at_ipart cart_propfar_at_ipart proc~cart_propfar_at_ipart->proc~fmm_l2l proc~fmm_solve_for_multipoles fmm_solve_for_multipoles proc~fmm_solve_for_multipoles->proc~tree_l2l proc~fmm_solve fmm_solve proc~fmm_solve->proc~tree_l2l proc~field_extd2d field_extD2D proc~field_extd2d->proc~cart_propfar_at_ipart proc~prepare_fmm_ext_ipd prepare_fmm_ext_ipd proc~field_extd2d->proc~prepare_fmm_ext_ipd proc~elec_prop_d2m elec_prop_D2M proc~elec_prop_d2m->proc~cart_propfar_at_ipart proc~prepare_fmm_ipd prepare_fmm_ipd proc~elec_prop_d2m->proc~prepare_fmm_ipd proc~elec_prop_m2m elec_prop_M2M proc~elec_prop_m2m->proc~cart_propfar_at_ipart proc~preapare_fmm_static preapare_fmm_static proc~elec_prop_m2m->proc~preapare_fmm_static proc~elec_prop_d2d elec_prop_D2D proc~elec_prop_d2d->proc~cart_propfar_at_ipart proc~elec_prop_d2d->proc~prepare_fmm_ipd proc~elec_prop_m2d elec_prop_M2D proc~elec_prop_m2d->proc~cart_propfar_at_ipart proc~elec_prop_m2d->proc~preapare_fmm_static proc~cart_prop_at_ipart cart_prop_at_ipart proc~cart_prop_at_ipart->proc~cart_propfar_at_ipart proc~prepare_fmm_ext_ipd->proc~fmm_solve_for_multipoles proc~preapare_fmm_static->proc~fmm_solve_for_multipoles proc~prepare_polelec prepare_polelec proc~prepare_polelec->proc~elec_prop_d2m proc~prepare_polelec->proc~elec_prop_d2d proc~prepare_polelec->proc~elec_prop_m2d proc~prepare_fixedelec prepare_fixedelec proc~prepare_fixedelec->proc~elec_prop_m2m proc~tmatvec_otf TMatVec_otf proc~tmatvec_otf->proc~field_extd2d proc~prepare_fmm_ipd->proc~prepare_fmm_ext_ipd proc~ommp_get_polelec_energy ommp_get_polelec_energy proc~ommp_get_polelec_energy->proc~prepare_polelec proc~energy_mm_pol energy_MM_pol proc~ommp_get_polelec_energy->proc~energy_mm_pol proc~ommp_set_external_field ommp_set_external_field proc~ommp_set_external_field->proc~prepare_polelec proc~polelec_geomgrad polelec_geomgrad proc~polelec_geomgrad->proc~prepare_polelec proc~energy_mm_pol->proc~prepare_polelec proc~fixedelec_geomgrad fixedelec_geomgrad proc~fixedelec_geomgrad->proc~prepare_fixedelec proc~energy_mm_mm energy_MM_MM proc~energy_mm_mm->proc~prepare_fixedelec proc~ommp_get_full_ele_energy ommp_get_full_ele_energy proc~ommp_get_full_ele_energy->proc~ommp_get_polelec_energy proc~ommp_get_fixedelec_energy ommp_get_fixedelec_energy proc~ommp_get_full_ele_energy->proc~ommp_get_fixedelec_energy proc~c_ommp_get_polelec_energy C_ommp_get_polelec_energy proc~c_ommp_get_polelec_energy->proc~ommp_get_polelec_energy proc~ommp_set_external_field_nomm ommp_set_external_field_nomm proc~ommp_set_external_field_nomm->proc~ommp_set_external_field proc~c_ommp_set_external_field C_ommp_set_external_field proc~c_ommp_set_external_field->proc~ommp_set_external_field proc~c_ommp_set_external_field_nomm C_ommp_set_external_field_nomm proc~c_ommp_set_external_field_nomm->proc~ommp_set_external_field proc~ommp_polelec_geomgrad ommp_polelec_geomgrad proc~ommp_polelec_geomgrad->proc~polelec_geomgrad proc~ommp_full_geomgrad ommp_full_geomgrad proc~ommp_full_geomgrad->proc~polelec_geomgrad proc~ommp_full_geomgrad->proc~fixedelec_geomgrad proc~ommp_fixedelec_geomgrad ommp_fixedelec_geomgrad proc~ommp_fixedelec_geomgrad->proc~fixedelec_geomgrad proc~ommp_get_fixedelec_energy->proc~energy_mm_mm proc~ommp_get_full_energy ommp_get_full_energy proc~ommp_get_full_energy->proc~ommp_get_full_ele_energy proc~c_ommp_polelec_geomgrad C_ommp_polelec_geomgrad proc~c_ommp_polelec_geomgrad->proc~ommp_polelec_geomgrad proc~c_ommp_full_geomgrad C_ommp_full_geomgrad proc~c_ommp_full_geomgrad->proc~ommp_full_geomgrad proc~c_ommp_fixedelec_geomgrad C_ommp_fixedelec_geomgrad proc~c_ommp_fixedelec_geomgrad->proc~ommp_fixedelec_geomgrad proc~c_ommp_get_full_ele_energy C_ommp_get_full_ele_energy proc~c_ommp_get_full_ele_energy->proc~ommp_get_full_ele_energy proc~c_ommp_get_fixedelec_energy C_ommp_get_fixedelec_energy proc~c_ommp_get_fixedelec_energy->proc~ommp_get_fixedelec_energy proc~c_ommp_get_full_energy C_ommp_get_full_energy proc~c_ommp_get_full_energy->proc~ommp_get_full_energy

Contents

Source Code

fmm_l2l_ztranslate_work

Source Code

subroutine fmm_l2l_ztranslate_work(z, p, vscales, vfact, alpha, &
    & src_l, beta, dst_l, work)
    ! Inputs
    integer, intent(in) :: p
    real(rp), intent(in) :: z, vscales((p+1)*(p+1)), &
        & vfact(2*p+1), alpha, src_l((p+1)*(p+1)), beta
    ! Output
    real(rp), intent(inout) :: dst_l((p+1)*(p+1))
    ! Temporary workspace
    real(rp), intent(out), target :: work(2*(p+1))
    ! Local variables
    real(rp) :: r1, r2, tmp1, tmp2
    integer :: j, k, n, indj, indn
    ! Pointers for temporary values of powers
    real(rp), pointer :: pow_r1(:), pow_r2(:)
    ! Init output
    if (abs(beta) < eps_rp) then
        dst_l = 0.0_rp
    else
        dst_l = beta * dst_l
    end if
    ! If harmonics have different centers
    if (abs(z) > eps_rp) then
        ! Prepare pointers
        pow_r1(1:p+1) => work(1:p+1)
        pow_r2(1:p+1) => work(p+2:2*(p+1))
        ! Get powers of r1 and r2
        r1 = z / 1.0_rp
        r2 = 1.0_rp / z
        pow_r1(1) = 1
        pow_r2(1) = 1
        do j = 2, p+1
            pow_r1(j) = pow_r1(j-1) * r1
            pow_r2(j) = pow_r2(j-1) * r2
        end do
        ! Do actual L2L
        do j = 0, p
            indj = j*j + j + 1
            do k = 0, j
                tmp1 = alpha * pow_r2(j+1) / vfact(j-k+1) / vfact(j+k+1) * &
                    & vscales(indj)
                do n = j, p
                    indn = n*n + n + 1
                    tmp2 = tmp1 * pow_r1(n+1) / vscales(indn) * &
                        & vfact(n-k+1) * vfact(n+k+1) / vfact(n-j+1) / &
                        & vfact(n-j+1)
                    if (mod(n+j, 2) .eq. 1) then
                        tmp2 = -tmp2
                    end if
                    if (k .eq. 0) then
                        dst_l(indj) = dst_l(indj) + tmp2*src_l(indn)
                    else
                        dst_l(indj+k) = dst_l(indj+k) + tmp2*src_l(indn+k)
                        dst_l(indj-k) = dst_l(indj-k) + tmp2*src_l(indn-k)
                    end if
                end do
            end do
        end do
    ! If harmonics are located at the same point
    else
        tmp1 = alpha
        do j = 0, p
            indj = j*j + j + 1
            do k = indj-j, indj+j
                dst_l(k) = dst_l(k) + src_l(k)*tmp1
            end do
            tmp1 = tmp1
        end do
    end if
end subroutine fmm_l2l_ztranslate_work