37subroutine ylmscale(p, vscales, v4pi2lp1, vscales_rel)
39 integer,
intent(in) :: p
41 real(dp),
intent(out) :: vscales((p+1)**2), v4pi2lp1(p+1), &
42 & vscales_rel((p+1)**2)
44 real(dp) :: tmp, twolp1
53 vscales_rel(ind) = tmp
54 vscales(ind) = one / tmp
56 tmp = vscales(ind) * sqrt2
59 tmp = -tmp / sqrt(dble((l-m+1)*(l+m)))
62 vscales_rel(ind+m) = tmp * v4pi2lp1(l+1)
63 vscales_rel(ind-m) = vscales_rel(ind+m)
81 & m2l_ztranslate_adj_coef)
83 integer,
intent(in) :: dmax, pm, pl
85 real(dp),
intent(out) :: vcnk((2*dmax+1)*(dmax+1)), &
86 & m2l_ztranslate_coef(pm+1, pl+1, pl+1), &
87 & m2l_ztranslate_adj_coef(pl+1, pl+1, pm+1)
89 integer :: i, indi, j, k, n, indjn
104 vcnk(indi+j) = vcnk(indi-i+j+1) + vcnk(indi-i+j)
112 vcnk(indi+j) = sqrt(vcnk(indi+j))
120 indjn = (j+n)*(j+n+1)/2 + 1
121 m2l_ztranslate_coef(n-k+1, k+1, j-k+1) = &
122 & tmp1 * vcnk(indjn+j-k) * vcnk(indjn+j+k)
123 m2l_ztranslate_adj_coef(j-k+1, k+1, n-k+1) = &
124 & m2l_ztranslate_coef(n-k+1, k+1, j-k+1)
144 real(dp),
intent(in) :: x(3)
146 real(dp),
intent(out) :: rho, ctheta, stheta, cphi, sphi
148 real(dp) :: max12, ssq12
150 if ((x(1) .eq. zero) .and. (x(2) .eq. zero))
then
152 ctheta = sign(one, x(3))
160 if (abs(x(2)) .gt. abs(x(1)))
then
162 ssq12 = one + (x(1)/x(2))**2
165 ssq12 = one + (x(2)/x(1))**2
167 stheta = max12 * sqrt(ssq12)
171 if (abs(x(3)) .gt. max12)
then
172 rho = one + ssq12*(max12/x(3))**2
173 rho = abs(x(3)) * sqrt(rho)
174 stheta = stheta / rho
177 rho = ssq12 + (x(3)/max12)**2
178 rho = max12 * sqrt(rho)
179 stheta = stheta / rho
214subroutine ylmbas(x, rho, ctheta, stheta, cphi, sphi, p, vscales, vylm, vplm, &
217 integer,
intent(in) :: p
218 real(dp),
intent(in) :: x(3)
219 real(dp),
intent(in) :: vscales((p+1)**2)
221 real(dp),
intent(out) :: rho, ctheta, stheta, cphi, sphi
222 real(dp),
intent(out) :: vylm((p+1)**2), vplm((p+1)**2)
223 real(dp),
intent(out) :: vcos(p+1), vsin(p+1)
226 real(dp) :: max12, ssq12, tmp
231 if (x(1) .eq. zero)
then
234 else if (abs(x(2)) .gt. abs(x(1)))
then
236 ssq12 = one + (x(1)/x(2))**2
239 ssq12 = one + (x(2)/x(1))**2
242 if (x(3) .eq. zero)
then
243 rho = max12 * sqrt(ssq12)
244 else if (abs(x(3)) .gt. max12)
then
245 rho = one + ssq12 *(max12/x(3))**2
246 rho = abs(x(3)) * sqrt(rho)
248 rho = ssq12 + (x(3)/max12)**2
249 rho = max12 * sqrt(rho)
252 if (rho .eq. zero)
then
256 stheta = max12 * sqrt(ssq12)
258 if (stheta .ne. zero)
then
264 stheta = stheta / rho
272 vylm(2) = - vscales(4) * stheta * sphi
273 vylm(3) = vscales(3) * ctheta
274 vylm(4) = - vscales(4) * stheta * cphi
281 call trgev(cphi, sphi, p, vcos, vsin)
286 vylm(2) = - vscales(4) * stheta * sphi
287 vylm(3) = vscales(3) * ctheta
288 vylm(4) = - vscales(4) * stheta * cphi
295 vylm(ind) = vscales(ind) * vplm(ind)
298 tmp = vplm(ind+m) * vscales(ind+m)
300 vylm(ind+m) = tmp * vcos(m+1)
302 vylm(ind-m) = tmp * vsin(m+1)
310 ctheta = sign(one, x(3))
324 vylm(ind) = vscales(ind)
330 vylm(ind) = ctheta * vscales(ind)
348subroutine trgev(cphi, sphi, p, vcos, vsin)
350 integer,
intent(in) :: p
351 real(dp),
intent(in) :: cphi, sphi
353 real(dp),
intent(out) :: vcos(p+1), vsin(p+1)
356 real(dp) :: c4phi, s4phi
380 vcos(3) = cphi**2 - sphi**2
381 vsin(3) = 2 * cphi * sphi
389 vcos(3) = cphi**2 - sphi**2
390 vsin(3) = 2 * cphi * sphi
391 vcos(4) = vcos(3)*cphi - vsin(3)*sphi
392 vsin(4) = vcos(3)*sphi + vsin(3)*cphi
400 vcos(3) = cphi**2 - sphi**2
401 vsin(3) = 2 * cphi * sphi
402 vcos(4) = vcos(3)*cphi - vsin(3)*sphi
403 vsin(4) = vcos(3)*sphi + vsin(3)*cphi
404 vcos(5) = vcos(3)**2 - vsin(3)**2
405 vsin(5) = 2 * vcos(3) * vsin(3)
411 vcos(m:m+3) = vcos(m-4:m-1)*c4phi - vsin(m-4:m-1)*s4phi
412 vsin(m:m+3) = vcos(m-4:m-1)*s4phi + vsin(m-4:m-1)*c4phi
417 vcos(p-1) = vcos(p-2)*cphi - vsin(p-2)*sphi
418 vsin(p-1) = vcos(p-2)*sphi + vsin(p-2)*cphi
419 vcos(p) = vcos(p-2)*vcos(3) - vsin(p-2)*vsin(3)
420 vsin(p) = vcos(p-2)*vsin(3) + vsin(p-2)*vcos(3)
421 vcos(p+1) = vcos(p-2)*vcos(4) - vsin(p-2)*vsin(4)
422 vsin(p+1) = vcos(p-2)*vsin(4) + vsin(p-2)*vcos(4)
424 vcos(p) = vcos(p-1)*cphi - vsin(p-1)*sphi
425 vsin(p) = vcos(p-1)*sphi + vsin(p-1)*cphi
426 vcos(p+1) = vcos(p-1)*vcos(3) - vsin(p-1)*vsin(3)
427 vsin(p+1) = vcos(p-1)*vsin(3) + vsin(p-1)*vcos(3)
429 vcos(p+1) = vcos(p)*cphi - vsin(p)*sphi
430 vsin(p+1) = vcos(p)*sphi + vsin(p)*cphi
452subroutine polleg(ctheta, stheta, p, vplm)
454 integer,
intent(in) :: p
455 real(dp),
intent(in) :: ctheta, stheta
457 real(dp),
intent(out) :: vplm((p+1)**2)
459 real(dp) :: work(p+1)
484 integer,
intent(in) :: p
485 real(dp),
intent(in) :: ctheta, stheta
487 real(dp),
intent(out) :: vplm((p+1)**2)
489 real(dp),
intent(out) :: work(p+1)
492 real(dp) :: tmp1, tmp2, tmp3, tmp4, pl2m, pl1m, plm, pmm
507 tmp1 = three * stheta
508 tmp2 = ctheta * ctheta
509 vplm(7) = 1.5d0*tmp2 - pt5
510 vplm(8) = -tmp1 * ctheta
511 vplm(9) = tmp1 * stheta
517 tmp1 = three * stheta
518 tmp2 = ctheta * ctheta
519 vplm(7) = 1.5d0*tmp2 - pt5
520 vplm(8) = -tmp1 * ctheta
521 vplm(9) = tmp1 * stheta
522 tmp3 = 2.5d0 * tmp2 - pt5
523 vplm(13) = tmp3*ctheta - ctheta
524 vplm(14) = -tmp1 * tmp3
526 vplm(15) = tmp4 * vplm(8)
527 vplm(16) = tmp4 * vplm(9)
533 work(l) = dble(2*l-1) * ctheta
544 tmp1 = dble(2*m+1) * pmm
550 vplm((m+1)*(m+3)) = pl1m
553 plm = work(l)*pl1m - dble(l+m-1)*pl2m
554 plm = plm / dble(l-m)
555 vplm(l*l+l+1+m) = plm
572subroutine modified_spherical_bessel_first_kind(lmax, argument, SI, DI, work)
575 integer,
intent(in) :: lmax
576 real(dp),
intent(in) :: argument
578 real(dp),
dimension(0:lmax),
intent(out) :: SI, DI
580 complex(dp),
dimension(0:max(1, lmax)),
intent(out) :: work
589 complex(dp) :: complex_argument
590 real(dp) :: scaling_factor
591 real(dp),
parameter :: fnu = 1.5d0
592 integer :: NZ, ierr, l
594 if (argument .eq. zero)
then
603 scaling_factor = sqrt(pi/(2*argument))
605 complex_argument = argument
607 call cbesi(complex_argument, pt5, 1, 1, work(0), nz, ierr)
608 if (ierr .ne. 0) stop
'Error in computing Bessel function of first kind'
610 if (lmax .gt. 0)
then
611 call cbesi(complex_argument, fnu, 1, lmax, work(1:lmax), nz, ierr)
612 if (ierr .ne. 0) stop
'Error in computing Bessel function of first kind'
615 call cbesi(complex_argument, fnu, 1, 1, work(1), nz, ierr)
616 if (ierr .ne. 0) stop
'Error in computing Bessel function of first kind'
619 si = real(work(0:lmax))
621 si = scaling_factor * si
623 di(0) = scaling_factor * real(work(1))
625 di(l)= si(l-1) - ((l+1.0d0)*si(l))/argument
627end subroutine modified_spherical_bessel_first_kind
635subroutine modified_spherical_bessel_second_kind(lmax, argument, SK, DK, work)
638 integer ,
intent(in) :: lmax
639 real(dp),
intent(in) :: argument
641 real(dp),
dimension(0:lmax),
intent(out) :: SK, DK
643 complex(dp),
dimension(0:max(1, lmax)),
intent(out) :: work
652 complex(dp) :: complex_argument
653 real(dp) :: scaling_factor
654 real(dp),
parameter :: fnu=1.5d0
655 integer :: NZ, ierr, l
658 scaling_factor = sqrt(2/(pi*argument))
660 complex_argument = argument
662 call cbesk(complex_argument, pt5, 1, 1, work(0), nz, ierr)
663 if (ierr .ne. 0) stop
'Error in computing Bessel function of second kind'
665 if (lmax .gt. 0)
then
666 call cbesk(complex_argument, fnu, 1, lmax, work(1:lmax), nz, ierr)
667 if (ierr .ne. 0) stop
'Error in computing Bessel function of second kind'
670 call cbesk(complex_argument, fnu, 1, 1, work(1), nz, ierr)
671 if (ierr .ne. 0) stop
'Error in computing Bessel function of second kind'
674 sk = real(work(0:lmax))
676 sk = scaling_factor * sk
678 dk(0) = -scaling_factor * real(work(1))
680 dk(l) = -sk(l-1) - ((l+1.0d0)*sk(l))/argument
682end subroutine modified_spherical_bessel_second_kind
710subroutine fmm_p2m(c, src_q, dst_r, p, vscales, beta, dst_m)
712 integer,
intent(in) :: p
713 real(dp),
intent(in) :: c(3), src_q, dst_r, vscales((p+1)**2), beta
715 real(dp),
intent(inout) :: dst_m((p+1)**2)
717 real(dp) :: work(2*(p+1)*(p+2))
719 call fmm_p2m_work(c, src_q, dst_r, p, vscales, beta, dst_m, work)
751 integer,
intent(in) :: p
752 real(dp),
intent(in) :: c(3), src_q, dst_r, vscales((p+1)**2), beta
754 real(dp),
intent(inout) :: dst_m((p+1)**2)
756 real(dp),
intent(out),
target :: work(2*(p+1)*(p+2))
758 real(dp) :: rho, ctheta, stheta, cphi, sphi, t, rcoef
759 integer :: n, ind, vylm1, vylm2, vplm1, vplm2, vcos1, vcos2, vsin1, vsin2
761 if (src_q .eq. zero)
then
763 if (beta .eq. zero)
then
783 call ylmbas(c, rho, ctheta, stheta, cphi, sphi, p, vscales, &
784 & work(vylm1:vylm2), work(vplm1:vplm2), work(vcos1:vcos2), &
787 if (rho .ne. zero)
then
791 if (beta .eq. zero)
then
795 dst_m(ind-n:ind+n) = t * work(ind-n:ind+n)
803 dst_m(ind-n:ind+n) = beta*dst_m(ind-n:ind+n) + &
804 & t*work(ind-n:ind+n)
811 if (beta .eq. zero)
then
812 dst_m(1) = src_q / dst_r / sqrt4pi
816 dst_m(1) = beta*dst_m(1) + src_q/dst_r/sqrt4pi
817 dst_m(2:) = beta * dst_m(2:)
848subroutine fmm_m2p(c, src_r, p, vscales_rel, alpha, src_m, beta, dst_v)
850 integer,
intent(in) :: p
851 real(dp),
intent(in) :: c(3), src_r, vscales_rel((p+1)*(p+1)), alpha, &
852 & src_m((p+1)*(p+1)), beta
854 real(dp),
intent(inout) :: dst_v
856 real(dp) :: work(p+1)
858 call fmm_m2p_work(c, src_r, p, vscales_rel, alpha, src_m, beta, dst_v, &
894 integer,
intent(in) :: p
895 real(dp),
intent(in) :: c(3), src_r, vscales_rel((p+1)*(p+1)), alpha, &
896 & src_m((p+1)*(p+1)), beta
898 real(dp),
intent(inout) :: dst_v
900 real(dp),
intent(out),
target :: work(p+1)
902 real(dp) :: rho, ctheta, stheta, cphi, sphi, rcoef, t, tmp, tmp1, tmp2, &
903 & tmp3, max12, ssq12, pl2m, pl1m, plm, pmm, cmphi, smphi, ylm
904 integer :: l, m, indl
906 if (beta .eq. zero)
then
913 if (alpha .eq. zero)
then
917 if (c(1) .eq. zero)
then
920 else if (abs(c(2)) .gt. abs(c(1)))
then
922 ssq12 = one + (c(1)/c(2))**2
925 ssq12 = one + (c(2)/c(1))**2
928 if (c(3) .eq. zero)
then
929 rho = max12 * sqrt(ssq12)
930 else if (abs(c(3)) .gt. max12)
then
931 rho = one + ssq12 *(max12/c(3))**2
932 rho = abs(c(3)) * sqrt(rho)
934 rho = ssq12 + (c(3)/max12)**2
935 rho = max12 * sqrt(rho)
939 if (rho .eq. zero)
then
945 stheta = max12 * sqrt(ssq12)
947 if (stheta .ne. zero)
then
953 stheta = stheta / rho
958 dst_v = dst_v + alpha*rcoef*vscales_rel(1)*src_m(1)
962 tmp = src_m(1) * vscales_rel(1)
964 tmp2 = ctheta * src_m(3) * vscales_rel(3)
965 tmp3 = vscales_rel(4) * stheta
966 tmp2 = tmp2 - tmp3*sphi*src_m(2)
967 tmp2 = tmp2 - tmp3*cphi*src_m(4)
968 dst_v = dst_v + alpha*rcoef*(tmp+rcoef*tmp2)
973 work(1) = alpha * rcoef
975 work(l+1) = rcoef * work(l)
980 dst_v = dst_v + work(1)*src_m(1)*vscales_rel(1)
985 ylm = pl1m * vscales_rel(3)
986 dst_v = dst_v + work(2)*src_m(3)*ylm
989 plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
991 ylm = plm * vscales_rel(l*l+l+1)
992 dst_v = dst_v + work(l+1)*src_m(l*l+l+1)*ylm
1004 ylm = pmm * vscales_rel((m+1)**2)
1005 tmp1 = cmphi*src_m((m+1)**2) + smphi*src_m(m*m+1)
1006 dst_v = dst_v + work(m+1)*ylm*tmp1
1008 tmp1 = dble(2*m+1) * pmm
1010 pmm = -stheta * tmp1
1012 pl1m = ctheta * tmp1
1013 ylm = pl1m * vscales_rel((m+1)*(m+3))
1014 tmp1 = cmphi*src_m((m+1)*(m+3)) + smphi*src_m((m+1)*(m+2)+1-m)
1015 dst_v = dst_v + work(m+2)*ylm*tmp1
1018 plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
1019 plm = plm / dble(l-m)
1020 ylm = plm * vscales_rel(l*l+l+1+m)
1021 tmp1 = cmphi*src_m(l*l+l+1+m) + smphi*src_m(l*l+l+1-m)
1022 dst_v = dst_v + work(l+1)*ylm*tmp1
1029 cmphi = cmphi*cphi - smphi*sphi
1030 smphi = tmp1*sphi + smphi*cphi
1033 ylm = pmm * vscales_rel((p+1)**2)
1034 tmp1 = cmphi*src_m((p+1)**2) + smphi*src_m(p*p+1)
1035 dst_v = dst_v + work(p+1)*ylm*tmp1
1044 if (c(3) .lt. zero)
then
1053 dst_v = dst_v + t*src_m(indl)*vscales_rel(indl)
1089 integer,
intent(in) :: p
1090 real(dp),
intent(in) :: c(3), src_r, vscales((p+1)*(p+1)), alpha, &
1091 & src_m((p+1)*(p+1)), beta
1093 real(dp),
intent(inout) :: dst_v
1095 complex(dp) :: work_complex(max(2, p+1))
1096 real(dp) :: work(p+1)
1098 real(dp) :: src_sk(p+1)
1100 call modified_spherical_bessel_second_kind(p, src_r, src_sk, work, &
1103 & dst_v, work_complex, work)
1133 & dst_v, work_complex, work)
1136 integer,
intent(in) :: p
1137 real(dp),
intent(in) :: c(3), vscales((p+1)*(p+1)), alpha, &
1138 & src_m((p+1)*(p+1)), beta, SK_ri(p+1)
1140 real(dp),
intent(inout) :: dst_v
1142 complex(dp),
intent(out) :: work_complex(p+1)
1143 real(dp),
intent(out) :: work(p+1)
1145 real(dp) :: rho, ctheta, stheta, cphi, sphi, tmp, tmp1, tmp2, &
1146 & tmp3, max12, ssq12, pl2m, pl1m, plm, pmm, cmphi, smphi, ylm, t
1147 real(dp) :: scaling_factor
1148 complex(dp) :: complex_argument
1149 integer :: l, m, indl, ierr, NZ
1150 real(dp),
parameter :: fnu=1.5d0
1152 if (beta .eq. zero)
then
1155 dst_v = beta * dst_v
1159 if (alpha .eq. zero)
then
1163 if (c(1) .eq. zero)
then
1166 else if (abs(c(2)) .gt. abs(c(1)))
then
1168 ssq12 = one + (c(1)/c(2))**2
1171 ssq12 = one + (c(2)/c(1))**2
1174 if (c(3) .eq. zero)
then
1175 rho = max12 * sqrt(ssq12)
1176 else if (abs(c(3)) .gt. max12)
then
1177 rho = one + ssq12 *(max12/c(3))**2
1178 rho = abs(c(3)) * sqrt(rho)
1180 rho = ssq12 + (c(3)/max12)**2
1181 rho = max12 * sqrt(rho)
1185 if (rho .eq. zero)
then
1189 scaling_factor = alpha * sqrt(two / (pi*rho))
1191 complex_argument = rho
1193 call cbesk(complex_argument, pt5, 1, 1, work_complex(1), nz, ierr)
1196 call cbesk(complex_argument, fnu, 1, p, work_complex(2:p+1), nz, ierr)
1198 work = scaling_factor * real(work_complex) / sk_ri
1201 stheta = max12 * sqrt(ssq12)
1203 if (stheta .ne. zero)
then
1205 cphi = c(1) / stheta
1206 sphi = c(2) / stheta
1209 stheta = stheta / rho
1214 dst_v = dst_v + work(1)*vscales(1)*src_m(1)
1218 tmp = work(1) * src_m(1) * vscales(1)
1220 tmp2 = ctheta * src_m(3) * vscales(3)
1221 tmp3 = stheta * vscales(4)
1222 tmp2 = tmp2 - tmp3*sphi*src_m(2)
1223 tmp2 = tmp2 - tmp3*cphi*src_m(4)
1224 dst_v = dst_v + tmp + work(2)*tmp2
1231 dst_v = dst_v + work(1)*src_m(1)*vscales(1)
1236 ylm = pl1m * vscales(3)
1237 dst_v = dst_v + work(2)*src_m(3)*ylm
1240 plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
1242 ylm = plm * vscales(l*l+l+1)
1243 dst_v = dst_v + work(l+1)*src_m(l*l+l+1)*ylm
1255 ylm = pmm * vscales((m+1)**2)
1256 tmp1 = cmphi*src_m((m+1)**2) + smphi*src_m(m*m+1)
1257 dst_v = dst_v + work(m+1)*ylm*tmp1
1259 tmp1 = dble(2*m+1) * pmm
1261 pmm = -stheta * tmp1
1263 pl1m = ctheta * tmp1
1264 ylm = pl1m * vscales((m+1)*(m+3))
1265 tmp1 = cmphi*src_m((m+1)*(m+3)) + smphi*src_m((m+1)*(m+2)+1-m)
1266 dst_v = dst_v + work(m+2)*ylm*tmp1
1269 plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
1270 plm = plm / dble(l-m)
1271 ylm = plm * vscales(l*l+l+1+m)
1272 tmp1 = cmphi*src_m(l*l+l+1+m) + smphi*src_m(l*l+l+1-m)
1273 dst_v = dst_v + work(l+1)*ylm*tmp1
1280 cmphi = cmphi*cphi - smphi*sphi
1281 smphi = tmp1*sphi + smphi*cphi
1284 ylm = pmm * vscales((p+1)**2)
1285 tmp1 = cmphi*src_m((p+1)**2) + smphi*src_m(p*p+1)
1286 dst_v = dst_v + work(p+1)*ylm*tmp1
1295 if (c(3) .lt. zero)
then
1302 dst_v = dst_v + t*work(l+1)*src_m(indl)*vscales(indl)
1312 dst_v = dst_v + work(l+1)*src_m(indl)*vscales(indl)
1318subroutine fmm_m2p_bessel_grad(c, src_r, p, vscales, alpha, src_m, beta, dst_g)
1321 integer,
intent(in) :: p
1322 real(dp),
intent(in) :: c(3), src_r, vscales((p+2)*(p+2)), alpha, &
1323 & src_m((p+1)*(p+1)), beta
1325 real(dp),
intent(inout) :: dst_g(3)
1327 complex(dp) :: work_complex(p+2)
1328 real(dp) :: work(p+2), src_sk(p+2), src_m_grad((p+2)**2, 3)
1330 call modified_spherical_bessel_second_kind(p+1, src_r, src_sk, work, &
1332 call fmm_m2m_bessel_grad(p, src_sk, vscales, src_m, src_m_grad)
1334 & beta, dst_g(1), work_complex, work)
1336 & beta, dst_g(2), work_complex, work)
1338 & beta, dst_g(3), work_complex, work)
1339end subroutine fmm_m2p_bessel_grad
1368 integer,
intent(in) :: p
1369 real(dp),
intent(in) :: c(3), src_q, dst_r, vscales_rel((p+1)*(p+1)), beta
1371 real(dp),
intent(inout) :: dst_m((p+1)**2)
1373 real(dp) :: work(p+1)
1406 integer,
intent(in) :: p
1407 real(dp),
intent(in) :: c(3), src_q, dst_r, vscales_rel((p+1)*(p+1)), beta
1409 real(dp),
intent(inout) :: dst_m((p+1)**2)
1411 real(dp),
intent(out),
target :: work(p+1)
1413 real(dp) :: rho, ctheta, stheta, cphi, sphi, rcoef, t, tmp, tmp1, tmp2, &
1414 & max12, ssq12, pl2m, pl1m, plm, pmm, cmphi, smphi, ylm
1415 integer :: l, m, indl
1417 if (src_q .eq. zero)
then
1419 if (beta .eq. zero)
then
1423 dst_m = beta * dst_m
1430 if (c(1) .eq. zero)
then
1433 else if (abs(c(2)) .gt. abs(c(1)))
then
1435 ssq12 = one + (c(1)/c(2))**2
1438 ssq12 = one + (c(2)/c(1))**2
1441 if (c(3) .eq. zero)
then
1442 rho = max12 * sqrt(ssq12)
1443 else if (abs(c(3)) .gt. max12)
then
1444 rho = one + ssq12 *(max12/c(3))**2
1445 rho = abs(c(3)) * sqrt(rho)
1447 rho = ssq12 + (c(3)/max12)**2
1448 rho = max12 * sqrt(rho)
1452 if (rho .eq. zero)
then
1458 stheta = max12 * sqrt(ssq12)
1460 if (stheta .ne. zero)
then
1462 cphi = c(1) / stheta
1463 sphi = c(2) / stheta
1466 stheta = stheta / rho
1471 if (beta .eq. zero)
then
1472 dst_m(1) = t * vscales_rel(1)
1474 dst_m(1) = beta*dst_m(1) + t*vscales_rel(1)
1478 if (beta .eq. zero)
then
1480 dst_m(1) = t * vscales_rel(1)
1483 tmp = -vscales_rel(4) * stheta
1485 dst_m(2) = tmp2 * sphi
1486 dst_m(3) = t * vscales_rel(3) * ctheta
1487 dst_m(4) = tmp2 * cphi
1490 dst_m(1) = beta*dst_m(1) + t*vscales_rel(1)
1493 tmp = -vscales_rel(4) * stheta
1495 dst_m(2) = beta*dst_m(2) + tmp2*sphi
1496 dst_m(3) = beta*dst_m(3) + t*vscales_rel(3)*ctheta
1497 dst_m(4) = beta*dst_m(4) + tmp2*cphi
1503 work(1) = src_q * rcoef
1505 work(l+1) = rcoef * work(l)
1508 if (beta .eq. zero)
then
1512 dst_m(1) = work(1) * vscales_rel(1)
1517 ylm = pl1m * vscales_rel(3)
1518 dst_m(3) = work(2) * ylm
1521 plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
1523 ylm = plm * vscales_rel(l*l+l+1)
1524 dst_m(l*l+l+1) = work(l+1) * ylm
1536 ylm = pmm * vscales_rel((m+1)**2)
1537 ylm = work(m+1) * ylm
1538 dst_m((m+1)**2) = cmphi * ylm
1539 dst_m(m*m+1) = smphi * ylm
1541 tmp1 = dble(2*m+1) * pmm
1543 pmm = -stheta * tmp1
1545 pl1m = ctheta * tmp1
1546 ylm = pl1m * vscales_rel((m+1)*(m+3))
1547 ylm = work(m+2) * ylm
1548 dst_m((m+1)*(m+3)) = cmphi * ylm
1549 dst_m((m+1)*(m+2)+1-m) = smphi * ylm
1552 plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
1553 plm = plm / dble(l-m)
1554 ylm = plm * vscales_rel(l*l+l+1+m)
1555 ylm = work(l+1) * ylm
1556 dst_m(l*l+l+1+m) = cmphi * ylm
1557 dst_m(l*l+l+1-m) = smphi * ylm
1564 cmphi = cmphi*cphi - smphi*sphi
1565 smphi = tmp1*sphi + smphi*cphi
1568 ylm = pmm * vscales_rel((p+1)**2)
1569 ylm = work(p+1) * ylm
1570 dst_m((p+1)**2) = cmphi * ylm
1571 dst_m(p*p+1) = smphi * ylm
1577 dst_m(1) = beta*dst_m(1) + work(1)*vscales_rel(1)
1582 ylm = pl1m * vscales_rel(3)
1583 dst_m(3) = beta*dst_m(3) + work(2)*ylm
1586 plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
1588 ylm = plm * vscales_rel(l*l+l+1)
1589 dst_m(l*l+l+1) = beta*dst_m(l*l+l+1) + work(l+1)*ylm
1601 ylm = pmm * vscales_rel((m+1)**2)
1602 ylm = work(m+1) * ylm
1603 dst_m((m+1)**2) = beta*dst_m((m+1)**2) + cmphi*ylm
1604 dst_m(m*m+1) = beta*dst_m(m*m+1) + smphi*ylm
1606 tmp1 = dble(2*m+1) * pmm
1608 pmm = -stheta * tmp1
1610 pl1m = ctheta * tmp1
1611 ylm = pl1m * vscales_rel((m+1)*(m+3))
1612 ylm = work(m+2) * ylm
1613 dst_m((m+1)*(m+3)) = beta*dst_m((m+1)*(m+3)) + cmphi*ylm
1614 dst_m((m+1)*(m+2)+1-m) = beta*dst_m((m+1)*(m+2)+1-m) + &
1618 plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
1619 plm = plm / dble(l-m)
1620 ylm = plm * vscales_rel(l*l+l+1+m)
1621 ylm = work(l+1) * ylm
1622 dst_m(l*l+l+1+m) = beta*dst_m(l*l+l+1+m) + cmphi*ylm
1623 dst_m(l*l+l+1-m) = beta*dst_m(l*l+l+1-m) + smphi*ylm
1630 cmphi = cmphi*cphi - smphi*sphi
1631 smphi = tmp1*sphi + smphi*cphi
1634 ylm = pmm * vscales_rel((p+1)**2)
1635 ylm = work(p+1) * ylm
1636 dst_m((p+1)**2) = beta*dst_m((p+1)**2) + cmphi*ylm
1637 dst_m(p*p+1) = beta*dst_m(p*p+1) + smphi*ylm
1647 if (c(3) .lt. zero)
then
1652 if (beta .eq. zero)
then
1657 dst_m(indl) = t * vscales_rel(indl)
1658 dst_m(indl-l:indl-1) = zero
1659 dst_m(indl+1:indl+l) = zero
1668 dst_m(indl) = beta*dst_m(indl) + t*vscales_rel(indl)
1669 dst_m(indl-l:indl-1) = beta * dst_m(indl-l:indl-1)
1670 dst_m(indl+1:indl+l) = beta * dst_m(indl+1:indl+l)
1705 integer,
intent(in) :: p
1706 real(dp),
intent(in) :: c(3), src_q, dst_r, vscales((p+1)*(p+1)), beta, &
1709 real(dp),
intent(inout) :: dst_m((p+1)**2)
1711 complex(dp) :: work_complex(max(2, p+1))
1712 real(dp) :: work(p+1)
1714 real(dp) :: dst_sk(p+1), ck(3)
1716 call modified_spherical_bessel_second_kind(p, kappa*dst_r, dst_sk, work, &
1720 & dst_m, work_complex, work)
1750 & work_complex, work)
1753 integer,
intent(in) :: p
1754 real(dp),
intent(in) :: c(3), src_q, dst_sk(p+1), vscales((p+1)*(p+1)), beta
1756 real(dp),
intent(inout) :: dst_m((p+1)**2)
1758 complex(dp),
intent(out) :: work_complex(p+1)
1759 real(dp),
intent(out) :: work(p+1)
1761 real(dp) :: rho, ctheta, stheta, cphi, sphi, t, tmp, tmp1, tmp2, &
1762 & max12, ssq12, pl2m, pl1m, plm, pmm, cmphi, smphi, ylm, &
1764 integer :: l, m, indl, NZ, ierr
1765 complex(dp) :: complex_argument
1767 if (src_q .eq. zero)
then
1769 if (beta .eq. zero)
then
1773 dst_m = beta * dst_m
1780 if (c(1) .eq. zero)
then
1783 else if (abs(c(2)) .gt. abs(c(1)))
then
1785 ssq12 = one + (c(1)/c(2))**2
1788 ssq12 = one + (c(2)/c(1))**2
1791 if (c(3) .eq. zero)
then
1792 rho = max12 * sqrt(ssq12)
1793 else if (abs(c(3)) .gt. max12)
then
1794 rho = one + ssq12 *(max12/c(3))**2
1795 rho = abs(c(3)) * sqrt(rho)
1797 rho = ssq12 + (c(3)/max12)**2
1798 rho = max12 * sqrt(rho)
1802 if (rho .eq. zero)
then
1806 scaling_factor = src_q * sqrt(two / (pi*rho))
1808 complex_argument = rho
1810 call cbesk(complex_argument, pt5, 1, 1, work_complex(1), nz, ierr)
1813 call cbesk(complex_argument, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
1815 work = scaling_factor * real(work_complex) / dst_sk
1817 stheta = max12 * sqrt(ssq12)
1819 if (stheta .ne. zero)
then
1821 cphi = c(1) / stheta
1822 sphi = c(2) / stheta
1825 stheta = stheta / rho
1829 if (beta .eq. zero)
then
1830 dst_m(1) = work(1) * vscales(1)
1832 dst_m(1) = beta*dst_m(1) + work(1)*vscales(1)
1836 if (beta .eq. zero)
then
1838 dst_m(1) = work(1) * vscales(1)
1840 tmp = -vscales(4) * stheta
1841 tmp2 = work(2) * tmp
1842 dst_m(2) = tmp2 * sphi
1843 dst_m(3) = work(2) * vscales(3) * ctheta
1844 dst_m(4) = tmp2 * cphi
1847 dst_m(1) = beta*dst_m(1) + work(1)*vscales(1)
1849 tmp = -vscales(4) * stheta
1850 tmp2 = work(2) * tmp
1851 dst_m(2) = beta*dst_m(2) + tmp2*sphi
1852 dst_m(3) = beta*dst_m(3) + work(2)*vscales(3)*ctheta
1853 dst_m(4) = beta*dst_m(4) + tmp2*cphi
1859 if (beta .eq. zero)
then
1863 dst_m(1) = work(1) * vscales(1)
1868 ylm = pl1m * vscales(3)
1869 dst_m(3) = work(2) * ylm
1872 plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
1874 ylm = plm * vscales(l*l+l+1)
1875 dst_m(l*l+l+1) = work(l+1) * ylm
1887 ylm = pmm * vscales((m+1)**2)
1888 ylm = work(m+1) * ylm
1889 dst_m((m+1)**2) = cmphi * ylm
1890 dst_m(m*m+1) = smphi * ylm
1892 tmp1 = dble(2*m+1) * pmm
1894 pmm = -stheta * tmp1
1896 pl1m = ctheta * tmp1
1897 ylm = pl1m * vscales((m+1)*(m+3))
1898 ylm = work(m+2) * ylm
1899 dst_m((m+1)*(m+3)) = cmphi * ylm
1900 dst_m((m+1)*(m+2)+1-m) = smphi * ylm
1903 plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
1904 plm = plm / dble(l-m)
1905 ylm = plm * vscales(l*l+l+1+m)
1906 ylm = work(l+1) * ylm
1907 dst_m(l*l+l+1+m) = cmphi * ylm
1908 dst_m(l*l+l+1-m) = smphi * ylm
1915 cmphi = cmphi*cphi - smphi*sphi
1916 smphi = tmp1*sphi + smphi*cphi
1919 ylm = pmm * vscales((p+1)**2)
1920 ylm = work(p+1) * ylm
1921 dst_m((p+1)**2) = cmphi * ylm
1922 dst_m(p*p+1) = smphi * ylm
1928 dst_m(1) = beta*dst_m(1) + work(1)*vscales(1)
1933 ylm = pl1m * vscales(3)
1934 dst_m(3) = beta*dst_m(3) + work(2)*ylm
1937 plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
1939 ylm = plm * vscales(l*l+l+1)
1940 dst_m(l*l+l+1) = beta*dst_m(l*l+l+1) + work(l+1)*ylm
1952 ylm = pmm * vscales((m+1)**2)
1953 ylm = work(m+1) * ylm
1954 dst_m((m+1)**2) = beta*dst_m((m+1)**2) + cmphi*ylm
1955 dst_m(m*m+1) = beta*dst_m(m*m+1) + smphi*ylm
1957 tmp1 = dble(2*m+1) * pmm
1959 pmm = -stheta * tmp1
1961 pl1m = ctheta * tmp1
1962 ylm = pl1m * vscales((m+1)*(m+3))
1963 ylm = work(m+2) * ylm
1964 dst_m((m+1)*(m+3)) = beta*dst_m((m+1)*(m+3)) + cmphi*ylm
1965 dst_m((m+1)*(m+2)+1-m) = beta*dst_m((m+1)*(m+2)+1-m) + &
1969 plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
1970 plm = plm / dble(l-m)
1971 ylm = plm * vscales(l*l+l+1+m)
1972 ylm = work(l+1) * ylm
1973 dst_m(l*l+l+1+m) = beta*dst_m(l*l+l+1+m) + cmphi*ylm
1974 dst_m(l*l+l+1-m) = beta*dst_m(l*l+l+1-m) + smphi*ylm
1981 cmphi = cmphi*cphi - smphi*sphi
1982 smphi = tmp1*sphi + smphi*cphi
1985 ylm = pmm * vscales((p+1)**2)
1986 ylm = work(p+1) * ylm
1987 dst_m((p+1)**2) = beta*dst_m((p+1)**2) + cmphi*ylm
1988 dst_m(p*p+1) = beta*dst_m(p*p+1) + smphi*ylm
1998 if (c(3) .lt. zero)
then
2001 if (beta .eq. zero)
then
2006 dst_m(indl) = t * work(l+1) * vscales(indl)
2007 dst_m(indl-l:indl-1) = zero
2008 dst_m(indl+1:indl+l) = zero
2017 dst_m(indl) = beta*dst_m(indl) + t*work(l+1)*vscales(indl)
2018 dst_m(indl-l:indl-1) = beta * dst_m(indl-l:indl-1)
2019 dst_m(indl+1:indl+l) = beta * dst_m(indl+1:indl+l)
2027 if (beta .eq. zero)
then
2032 dst_m(indl) = work(l+1) * vscales(indl)
2033 dst_m(indl-l:indl-1) = zero
2034 dst_m(indl+1:indl+l) = zero
2041 dst_m(indl) = beta*dst_m(indl) + work(l+1)*vscales(indl)
2042 dst_m(indl-l:indl-1) = beta * dst_m(indl-l:indl-1)
2043 dst_m(indl+1:indl+l) = beta * dst_m(indl+1:indl+l)
2075 integer,
intent(in) :: p
2076 real(dp),
intent(in) :: c(3), r, vscales((p+1)**2)
2078 real(dp),
intent(out) :: mat((p+1)**2)
2080 real(dp) :: rho, ctheta, stheta, cphi, sphi, vcos(p+1), vsin(p+1)
2081 real(dp) :: vylm((p+1)**2), vplm((p+1)**2), rcoef, t
2084 call ylmbas(c, rho, ctheta, stheta, cphi, sphi, p, vscales, vylm, vplm, &
2088 if (rho .eq. zero)
then
2097 mat(ind-n:ind+n) = t / vscales(ind)**2 * vylm(ind-n:ind+n)
2126subroutine fmm_l2p(c, src_r, p, vscales, alpha, src_l, beta, dst_v)
2128 integer,
intent(in) :: p
2129 real(dp),
intent(in) :: c(3), src_r, vscales((p+1)*(p+1)), alpha, &
2130 & src_l((p+1)*(p+1)), beta
2132 real(dp),
intent(inout) :: dst_v
2134 real(dp) :: work(p+1)
2136 call fmm_l2p_work(c, src_r, p, vscales, alpha, src_l, beta, dst_v, work)
2165subroutine fmm_l2p_work(c, src_r, p, vscales_rel, alpha, src_l, beta, dst_v, &
2168 integer,
intent(in) :: p
2169 real(dp),
intent(in) :: c(3), src_r, vscales_rel((p+1)*(p+1)), alpha, &
2170 & src_l((p+1)*(p+1)), beta
2172 real(dp),
intent(inout) :: dst_v
2174 real(dp),
intent(out),
target :: work(p+1)
2176 real(dp) :: rho, ctheta, stheta, cphi, sphi, rcoef, t, tmp, tmp1, tmp2, &
2177 & tmp3, max12, ssq12, pl2m, pl1m, plm, pmm, cmphi, smphi, ylm
2178 integer :: l, m, indl
2180 if (beta .eq. zero)
then
2183 dst_v = beta * dst_v
2187 if (alpha .eq. zero)
then
2191 if (c(1) .eq. zero)
then
2194 else if (abs(c(2)) .gt. abs(c(1)))
then
2196 ssq12 = one + (c(1)/c(2))**2
2199 ssq12 = one + (c(2)/c(1))**2
2202 if (c(3) .eq. zero)
then
2203 rho = max12 * sqrt(ssq12)
2204 else if (abs(c(3)) .gt. max12)
then
2205 rho = one + ssq12 *(max12/c(3))**2
2206 rho = abs(c(3)) * sqrt(rho)
2208 rho = ssq12 + (c(3)/max12)**2
2209 rho = max12 * sqrt(rho)
2213 if (rho .eq. zero)
then
2214 dst_v = dst_v + alpha*src_l(1)*vscales_rel(1)
2220 stheta = max12 * sqrt(ssq12)
2222 if (stheta .ne. zero)
then
2224 cphi = c(1) / stheta
2225 sphi = c(2) / stheta
2228 stheta = stheta / rho
2233 dst_v = dst_v + alpha*vscales_rel(1)*src_l(1)
2237 tmp = src_l(1) * vscales_rel(1)
2239 tmp2 = ctheta * src_l(3) * vscales_rel(3)
2240 tmp3 = vscales_rel(4) * stheta
2241 tmp2 = tmp2 - tmp3*sphi*src_l(2)
2242 tmp2 = tmp2 - tmp3*cphi*src_l(4)
2243 dst_v = dst_v + alpha*(tmp+rcoef*tmp2)
2250 work(l+1) = rcoef * work(l)
2255 dst_v = dst_v + work(1)*src_l(1)*vscales_rel(1)
2260 ylm = pl1m * vscales_rel(3)
2261 dst_v = dst_v + work(2)*src_l(3)*ylm
2264 plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
2266 ylm = plm * vscales_rel(l*l+l+1)
2267 dst_v = dst_v + work(l+1)*src_l(l*l+l+1)*ylm
2279 ylm = pmm * vscales_rel((m+1)**2)
2280 tmp1 = cmphi*src_l((m+1)**2) + smphi*src_l(m*m+1)
2281 dst_v = dst_v + work(m+1)*ylm*tmp1
2283 tmp1 = dble(2*m+1) * pmm
2285 pmm = -stheta * tmp1
2287 pl1m = ctheta * tmp1
2288 ylm = pl1m * vscales_rel((m+1)*(m+3))
2289 tmp1 = cmphi*src_l((m+1)*(m+3)) + smphi*src_l((m+1)*(m+2)+1-m)
2290 dst_v = dst_v + work(m+2)*ylm*tmp1
2293 plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
2294 plm = plm / dble(l-m)
2295 ylm = plm * vscales_rel(l*l+l+1+m)
2296 tmp1 = cmphi*src_l(l*l+l+1+m) + smphi*src_l(l*l+l+1-m)
2297 dst_v = dst_v + work(l+1)*ylm*tmp1
2304 cmphi = cmphi*cphi - smphi*sphi
2305 smphi = tmp1*sphi + smphi*cphi
2308 ylm = pmm * vscales_rel((p+1)**2)
2309 tmp1 = cmphi*src_l((p+1)**2) + smphi*src_l(p*p+1)
2310 dst_v = dst_v + work(p+1)*ylm*tmp1
2316 rcoef = sign(rcoef, c(3))
2324 dst_v = dst_v + t*src_l(indl)*vscales_rel(indl)
2358 integer,
intent(in) :: p
2359 real(dp),
intent(in) :: c(3), src_r, vscales((p+1)*(p+1)), alpha, &
2360 & src_l((p+1)*(p+1)), beta
2362 real(dp),
intent(inout) :: dst_v
2364 complex(dp) :: work_complex(max(2, p+1))
2365 real(dp) :: work(p+1)
2367 real(dp) :: src_si(p+1)
2369 call modified_spherical_bessel_first_kind(p, src_r, src_si, work, &
2372 & dst_v, work_complex, work)
2402 & dst_v, work_complex, work)
2405 integer,
intent(in) :: p
2406 real(dp),
intent(in) :: c(3), vscales((p+1)*(p+1)), alpha, &
2407 & src_l((p+1)*(p+1)), beta, SI_ri(p+1)
2409 real(dp),
intent(inout) :: dst_v
2411 complex(dp),
intent(out) :: work_complex(p+1)
2412 real(dp),
intent(out) :: work(p+1)
2414 real(dp) :: rho, ctheta, stheta, cphi, sphi, t, tmp, tmp1, tmp2, &
2415 & tmp3, max12, ssq12, pl2m, pl1m, plm, pmm, cmphi, smphi, ylm
2416 real(dp) :: scaling_factor
2417 complex(dp) :: complex_argument
2418 integer :: l, m, indl, ierr, NZ
2419 real(dp),
parameter :: fnu=1.5d0
2421 if (beta .eq. zero)
then
2424 dst_v = beta * dst_v
2428 if (alpha .eq. zero)
then
2432 if (c(1) .eq. zero)
then
2435 else if (abs(c(2)) .gt. abs(c(1)))
then
2437 ssq12 = one + (c(1)/c(2))**2
2440 ssq12 = one + (c(2)/c(1))**2
2443 if (c(3) .eq. zero)
then
2444 rho = max12 * sqrt(ssq12)
2445 else if (abs(c(3)) .gt. max12)
then
2446 rho = one + ssq12 *(max12/c(3))**2
2447 rho = abs(c(3)) * sqrt(rho)
2449 rho = ssq12 + (c(3)/max12)**2
2450 rho = max12 * sqrt(rho)
2454 if (rho .eq. zero)
then
2455 dst_v = dst_v + alpha*src_l(1)*vscales(1)/si_ri(1)
2459 scaling_factor = alpha * sqrt(pi / (2*rho))
2461 complex_argument = rho
2463 call cbesi(complex_argument, pt5, 1, 1, work_complex(1), nz, ierr)
2466 call cbesi(complex_argument, fnu, 1, p, work_complex(2:p+1), nz, ierr)
2468 work = scaling_factor * real(work_complex) / si_ri
2471 stheta = max12 * sqrt(ssq12)
2473 if (stheta .ne. zero)
then
2475 cphi = c(1) / stheta
2476 sphi = c(2) / stheta
2479 stheta = stheta / rho
2484 dst_v = dst_v + work(1)*vscales(1)*src_l(1)
2488 tmp = work(1) * src_l(1) * vscales(1)
2490 tmp2 = ctheta * src_l(3) * vscales(3)
2491 tmp3 = stheta * vscales(4)
2492 tmp2 = tmp2 - tmp3*sphi*src_l(2)
2493 tmp2 = tmp2 - tmp3*cphi*src_l(4)
2494 dst_v = dst_v + tmp + work(2)*tmp2
2501 dst_v = dst_v + work(1)*src_l(1)*vscales(1)
2506 ylm = pl1m * vscales(3)
2507 dst_v = dst_v + work(2)*src_l(3)*ylm
2510 plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
2512 ylm = plm * vscales(l*l+l+1)
2513 dst_v = dst_v + work(l+1)*src_l(l*l+l+1)*ylm
2525 ylm = pmm * vscales((m+1)**2)
2526 tmp1 = cmphi*src_l((m+1)**2) + smphi*src_l(m*m+1)
2527 dst_v = dst_v + work(m+1)*ylm*tmp1
2529 tmp1 = dble(2*m+1) * pmm
2531 pmm = -stheta * tmp1
2533 pl1m = ctheta * tmp1
2534 ylm = pl1m * vscales((m+1)*(m+3))
2535 tmp1 = cmphi*src_l((m+1)*(m+3)) + smphi*src_l((m+1)*(m+2)+1-m)
2536 dst_v = dst_v + work(m+2)*ylm*tmp1
2539 plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
2540 plm = plm / dble(l-m)
2541 ylm = plm * vscales(l*l+l+1+m)
2542 tmp1 = cmphi*src_l(l*l+l+1+m) + smphi*src_l(l*l+l+1-m)
2543 dst_v = dst_v + work(l+1)*ylm*tmp1
2550 cmphi = cmphi*cphi - smphi*sphi
2551 smphi = tmp1*sphi + smphi*cphi
2554 ylm = pmm * vscales((p+1)**2)
2555 tmp1 = cmphi*src_l((p+1)**2) + smphi*src_l(p*p+1)
2556 dst_v = dst_v + work(p+1)*ylm*tmp1
2565 if (c(3) .lt. zero)
then
2572 dst_v = dst_v + t*work(l+1)*src_l(indl)*vscales(indl)
2582 dst_v = dst_v + work(l+1)*src_l(indl)*vscales(indl)
2588subroutine fmm_l2p_bessel_grad(c, src_r, p, vscales, alpha, src_l, beta, dst_g)
2591 integer,
intent(in) :: p
2592 real(dp),
intent(in) :: c(3), src_r, vscales((p+2)*(p+2)), alpha, &
2593 & src_l((p+1)*(p+1)), beta
2595 real(dp),
intent(inout) :: dst_g(3)
2597 complex(dp) :: work_complex(p+2)
2598 real(dp) :: work(p+2), src_si(p+2), src_l_grad((p+2)**2, 3)
2600 call modified_spherical_bessel_first_kind(p+1, src_r, src_si, work, &
2602 call fmm_l2l_bessel_grad(p, src_si, vscales, src_l, src_l_grad)
2604 & beta, dst_g(1), work_complex, work)
2606 & beta, dst_g(2), work_complex, work)
2608 & beta, dst_g(3), work_complex, work)
2609end subroutine fmm_l2p_bessel_grad
2637 integer,
intent(in) :: p
2638 real(dp),
intent(in) :: c(3), src_q, dst_r, vscales_rel((p+1)*(p+1)), beta
2640 real(dp),
intent(inout) :: dst_l((p+1)**2)
2642 real(dp) :: work(p+1)
2674 integer,
intent(in) :: p
2675 real(dp),
intent(in) :: c(3), src_q, dst_r, vscales_rel((p+1)*(p+1)), beta
2677 real(dp),
intent(inout) :: dst_l((p+1)**2)
2679 real(dp),
intent(out),
target :: work(p+1)
2681 real(dp) :: rho, ctheta, stheta, cphi, sphi, rcoef, t, tmp, tmp1, tmp2, &
2682 & max12, ssq12, pl2m, pl1m, plm, pmm, cmphi, smphi, ylm
2683 integer :: l, m, indl
2685 if (src_q .eq. zero)
then
2687 if (beta .eq. zero)
then
2691 dst_l = beta * dst_l
2698 if (c(1) .eq. zero)
then
2701 else if (abs(c(2)) .gt. abs(c(1)))
then
2703 ssq12 = one + (c(1)/c(2))**2
2706 ssq12 = one + (c(2)/c(1))**2
2709 if (c(3) .eq. zero)
then
2710 rho = max12 * sqrt(ssq12)
2711 else if (abs(c(3)) .gt. max12)
then
2712 rho = one + ssq12 *(max12/c(3))**2
2713 rho = abs(c(3)) * sqrt(rho)
2715 rho = ssq12 + (c(3)/max12)**2
2716 rho = max12 * sqrt(rho)
2720 if (rho .eq. zero)
then
2721 if (beta .eq. zero)
then
2722 dst_l(1) = src_q * vscales_rel(1)
2725 dst_l(1) = beta*dst_l(1) + src_q*vscales_rel(1)
2726 dst_l(2:) = beta * dst_l(2:)
2733 stheta = max12 * sqrt(ssq12)
2735 if (stheta .ne. zero)
then
2737 cphi = c(1) / stheta
2738 sphi = c(2) / stheta
2741 stheta = stheta / rho
2746 if (beta .eq. zero)
then
2747 dst_l(1) = src_q * vscales_rel(1)
2749 dst_l(1) = beta*dst_l(1) + src_q*vscales_rel(1)
2755 tmp2 = -tmp * vscales_rel(4) * stheta
2756 if (beta .eq. zero)
then
2757 dst_l(1) = src_q * vscales_rel(1)
2758 dst_l(2) = tmp2 * sphi
2759 dst_l(3) = tmp * vscales_rel(3) * ctheta
2760 dst_l(4) = tmp2 * cphi
2762 dst_l(1) = beta*dst_l(1) + src_q*vscales_rel(1)
2763 dst_l(2) = beta*dst_l(2) + tmp2*sphi
2764 dst_l(3) = beta*dst_l(3) + tmp*vscales_rel(3)*ctheta
2765 dst_l(4) = beta*dst_l(4) + tmp2*cphi
2773 work(l+1) = rcoef * work(l)
2776 if (beta .eq. zero)
then
2780 dst_l(1) = work(1) * vscales_rel(1)
2785 ylm = pl1m * vscales_rel(3)
2786 dst_l(3) = work(2) * ylm
2789 plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
2791 ylm = plm * vscales_rel(l*l+l+1)
2792 dst_l(l*l+l+1) = work(l+1) * ylm
2804 ylm = pmm * vscales_rel((m+1)**2)
2805 ylm = work(m+1) * ylm
2806 dst_l((m+1)**2) = cmphi * ylm
2807 dst_l(m*m+1) = smphi * ylm
2809 tmp1 = dble(2*m+1) * pmm
2811 pmm = -stheta * tmp1
2813 pl1m = ctheta * tmp1
2814 ylm = pl1m * vscales_rel((m+1)*(m+3))
2815 ylm = work(m+2) * ylm
2816 dst_l((m+1)*(m+3)) = cmphi * ylm
2817 dst_l((m+1)*(m+2)+1-m) = smphi * ylm
2820 plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
2821 plm = plm / dble(l-m)
2822 ylm = plm * vscales_rel(l*l+l+1+m)
2823 ylm = work(l+1) * ylm
2824 dst_l(l*l+l+1+m) = cmphi * ylm
2825 dst_l(l*l+l+1-m) = smphi * ylm
2832 cmphi = cmphi*cphi - smphi*sphi
2833 smphi = tmp1*sphi + smphi*cphi
2836 ylm = pmm * vscales_rel((p+1)**2)
2837 ylm = work(p+1) * ylm
2838 dst_l((p+1)**2) = cmphi * ylm
2839 dst_l(p*p+1) = smphi * ylm
2845 dst_l(1) = beta*dst_l(1) + work(1)*vscales_rel(1)
2850 ylm = pl1m * vscales_rel(3)
2851 dst_l(3) = beta*dst_l(3) + work(2)*ylm
2854 plm = dble(2*l-1)*ctheta*pl1m - dble(l-1)*pl2m
2856 ylm = plm * vscales_rel(l*l+l+1)
2857 dst_l(l*l+l+1) = beta*dst_l(l*l+l+1) + work(l+1)*ylm
2869 ylm = pmm * vscales_rel((m+1)**2)
2870 ylm = work(m+1) * ylm
2871 dst_l((m+1)**2) = beta*dst_l((m+1)**2) + cmphi*ylm
2872 dst_l(m*m+1) = beta*dst_l(m*m+1) + smphi*ylm
2874 tmp1 = dble(2*m+1) * pmm
2876 pmm = -stheta * tmp1
2878 pl1m = ctheta * tmp1
2879 ylm = pl1m * vscales_rel((m+1)*(m+3))
2880 ylm = work(m+2) * ylm
2881 dst_l((m+1)*(m+3)) = beta*dst_l((m+1)*(m+3)) + cmphi*ylm
2882 dst_l((m+1)*(m+2)+1-m) = beta*dst_l((m+1)*(m+2)+1-m) + &
2886 plm = dble(2*l-1)*ctheta*pl1m - dble(l+m-1)*pl2m
2887 plm = plm / dble(l-m)
2888 ylm = plm * vscales_rel(l*l+l+1+m)
2889 ylm = work(l+1) * ylm
2890 dst_l(l*l+l+1+m) = beta*dst_l(l*l+l+1+m) + cmphi*ylm
2891 dst_l(l*l+l+1-m) = beta*dst_l(l*l+l+1-m) + smphi*ylm
2898 cmphi = cmphi*cphi - smphi*sphi
2899 smphi = tmp1*sphi + smphi*cphi
2902 ylm = pmm * vscales_rel((p+1)**2)
2903 ylm = work(p+1) * ylm
2904 dst_l((p+1)**2) = beta*dst_l((p+1)**2) + cmphi*ylm
2905 dst_l(p*p+1) = beta*dst_l(p*p+1) + smphi*ylm
2912 rcoef = sign(rcoef, c(3))
2917 if (beta .eq. zero)
then
2922 dst_l(indl) = t * vscales_rel(indl)
2923 dst_l(indl-l:indl-1) = zero
2924 dst_l(indl+1:indl+l) = zero
2934 dst_l(indl) = beta*dst_l(indl) + t*vscales_rel(indl)
2935 dst_l(indl-l:indl-1) = beta * dst_l(indl-l:indl-1)
2936 dst_l(indl+1:indl+l) = beta * dst_l(indl+1:indl+l)
3000 integer,
intent(in) :: p
3001 real(dp),
intent(in) :: r1(-1:1, -1:1), alpha, src((p+1)*(p+1)), beta
3003 real(dp),
intent(inout) :: dst((p+1)*(p+1))
3005 real(dp) :: work(2*(2*p+1)*(2*p+3))
3025 integer,
intent(in) :: p
3026 real(dp),
intent(in) :: r1(-1:1, -1:1), alpha, src((p+1)*(p+1)), beta
3028 real(dp),
intent(out) :: dst((p+1)*(p+1))
3030 real(dp),
intent(out),
target :: work(2*(2*p+1)*(2*p+3))
3033 integer :: l, m, n, ind
3035 real(dp),
pointer :: r_prev(:, :), r(:, :), scal_uvw_m(:), scal_u_n(:), &
3036 & scal_v_n(:), scal_w_n(:), r_swap(:, :)
3038 if (alpha .eq. zero)
then
3040 if (beta .eq. zero)
then
3050 if (beta .eq. zero)
then
3053 dst(1) = alpha * src(1)
3058 dst(2) = alpha * dot_product(r1(-1:1,-1), src(2:4))
3059 dst(3) = alpha * dot_product(r1(-1:1,0), src(2:4))
3060 dst(4) = alpha * dot_product(r1(-1:1,1), src(2:4))
3067 r_prev(-p:p, -p:p) => work(1:l)
3069 r(-p:p, -p:p) => work(l+1:m)
3071 scal_uvw_m(-p:p) => work(m+1:l)
3073 scal_u_n(-p:p) => work(l+1:m)
3075 scal_v_n(-p:p) => work(m+1:l)
3077 scal_w_n(-p:p) => work(l+1:m)
3079 r(2, 2) = (r1(1, 1)*r1(1, 1) - r1(1, -1)*r1(1, -1) - &
3080 & r1(-1, 1)*r1(-1, 1) + r1(-1, -1)*r1(-1, -1)) / two
3081 r(2, 1) = r1(1, 1)*r1(1, 0) - r1(-1, 1)*r1(-1, 0)
3082 r(2, 0) = sqrt3 / two * (r1(1, 0)*r1(1, 0) - r1(-1, 0)*r1(-1, 0))
3083 r(2, -1) = r1(1, -1)*r1(1, 0) - r1(-1, -1)*r1(-1, 0)
3084 r(2, -2) = r1(1, 1)*r1(1, -1) - r1(-1, 1)*r1(-1, -1)
3085 r(1, 2) = r1(1, 1)*r1(0, 1) - r1(1, -1)*r1(0, -1)
3086 r(1, 1) = r1(1, 1)*r1(0, 0) + r1(1, 0)*r1(0, 1)
3087 r(1, 0) = sqrt3 * r1(1, 0) * r1(0, 0)
3088 r(1, -1) = r1(1, -1)*r1(0, 0) + r1(1, 0)*r1(0, -1)
3089 r(1, -2) = r1(1, 1)*r1(0, -1) + r1(1, -1)*r1(0, 1)
3090 r(0, 2) = sqrt3 / two * (r1(0, 1)*r1(0, 1) - r1(0, -1)*r1(0, -1))
3091 r(0, 1) = sqrt3 * r1(0, 1) * r1(0, 0)
3092 r(0, 0) = (three*r1(0, 0)*r1(0, 0)-one) / two
3093 r(0, -1) = sqrt3 * r1(0, -1) * r1(0, 0)
3094 r(0, -2) = sqrt3 * r1(0, 1) * r1(0, -1)
3095 r(-1, 2) = r1(-1, 1)*r1(0, 1) - r1(-1, -1)*r1(0, -1)
3096 r(-1, 1) = r1(-1, 1)*r1(0, 0) + r1(-1, 0)*r1(0, 1)
3097 r(-1, 0) = sqrt3 * r1(-1, 0) * r1(0, 0)
3098 r(-1, -1) = r1(-1, -1)*r1(0, 0) + r1(0, -1)*r1(-1, 0)
3099 r(-1, -2) = r1(-1, 1)*r1(0, -1) + r1(-1, -1)*r1(0, 1)
3100 r(-2, 2) = r1(1, 1)*r1(-1, 1) - r1(1, -1)*r1(-1, -1)
3101 r(-2, 1) = r1(1, 1)*r1(-1, 0) + r1(1, 0)*r1(-1, 1)
3102 r(-2, 0) = sqrt3 * r1(1, 0) * r1(-1, 0)
3103 r(-2, -1) = r1(1, -1)*r1(-1, 0) + r1(1, 0)*r1(-1, -1)
3104 r(-2, -2) = r1(1, 1)*r1(-1, -1) + r1(1, -1)*r1(-1, 1)
3106 dst(7+m) = alpha * dot_product(r(-2:2, m), src(5:9))
3115 scal_uvw_m(0) = dble(l)
3117 u = sqrt(dble(l*l-m*m))
3122 u = sqrt(dble(u*(u-one)))
3125 scal_u_n(0) = dble(l)
3126 scal_v_n(0) = -sqrt(dble(l*(l-1))) / sqrt2
3129 u = sqrt(dble(l*l-n*n))
3133 v = sqrt(v*(v-one)) / two
3137 w = -sqrt(w*(w-one)) / two
3141 u = sqrt(dble(2*l-1))
3146 v = sqrt(dble((2*l-1)*(2*l-2))) / two
3149 v = sqrt(dble(2*l*(2*l-1))) / two
3152 scal_w_n(l-1) = zero
3155 scal_w_n(1-l) = zero
3158 v = r1(1, 1)*r_prev(l-1, l-1) - r1(1, -1)*r_prev(l-1, 1-l) - &
3159 & r1(-1, 1)*r_prev(1-l, l-1) + r1(-1, -1)*r_prev(1-l, 1-l)
3160 r(l, l) = v * scal_v_n(l) / scal_uvw_m(l)
3163 v = r1(1, 1)*r_prev(1-l, l-1) - r1(1, -1)*r_prev(1-l, 1-l) + &
3164 & r1(-1, 1)*r_prev(l-1, l-1) - r1(-1, -1)*r_prev(l-1, 1-l)
3165 r(-l, l) = v * scal_v_n(l) / scal_uvw_m(l)
3168 u = r1(0, 1)*r_prev(l-1, l-1) - r1(0, -1)*r_prev(l-1, 1-l)
3169 v = r1(1, 1)*r_prev(l-2, l-1) - r1(1, -1)*r_prev(l-2, 1-l) - &
3170 & r1(-1, 1)*r_prev(2-l, l-1) + r1(-1, -1)*r_prev(2-l, 1-l)
3171 r(l-1, l) = (u*scal_u_n(l-1)+v*scal_v_n(l-1)) / scal_uvw_m(l)
3174 u = r1(0, 1)*r_prev(1-l, l-1) - r1(0, -1)*r_prev(1-l, 1-l)
3175 v = r1(1, 1)*r_prev(2-l, l-1) - r1(1, -1)*r_prev(2-l, 1-l) + &
3176 & r1(-1, 1)*r_prev(l-2, l-1) - r1(-1, -1)*r_prev(l-2, 1-l)
3177 r(1-l, l) = (u*scal_u_n(l-1)+v*scal_v_n(l-1)) / scal_uvw_m(l)
3180 u = r1(0, 1)*r_prev(1, l-1) - r1(0, -1)*r_prev(1, 1-l)
3181 v = r1(1, 1)*r_prev(0, l-1) - r1(1, -1)*r_prev(0, 1-l)
3182 w = r1(1, 1)*r_prev(2, l-1) - r1(1, -1)*r_prev(2, 1-l) + &
3183 & r1(-1, 1)*r_prev(-2, l-1) - r1(-1, -1)*r_prev(-2, 1-l)
3184 r(1, l) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
3185 r(1, l) = r(1, l) / scal_uvw_m(l)
3188 u = r1(0, 1)*r_prev(-1, l-1) - r1(0, -1)*r_prev(-1, 1-l)
3189 v = r1(-1, 1)*r_prev(0, l-1) - r1(-1, -1)*r_prev(0, 1-l)
3190 w = r1(1, 1)*r_prev(-2, l-1) - r1(1, -1)*r_prev(-2, 1-l) - &
3191 & r1(-1, 1)*r_prev(2, l-1) + r1(-1, -1)*r_prev(2, 1-l)
3192 r(-1, l) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
3193 r(-1, l) = r(-1, l) / scal_uvw_m(l)
3196 u = r1(0, 1)*r_prev(0, l-1) - r1(0, -1)*r_prev(0, 1-l)
3197 v = r1(1, 1)*r_prev(1, l-1) - r1(1, -1)*r_prev(1, 1-l) + &
3198 & r1(-1, 1)*r_prev(-1, l-1) - r1(-1, -1)*r_prev(-1, 1-l)
3199 r(0, l) = (u*scal_u_n(0) + v*scal_v_n(0)) / scal_uvw_m(l)
3203 u = r1(0, 1)*r_prev(n, l-1) - r1(0, -1)*r_prev(n, 1-l)
3204 v = r1(1, 1)*r_prev(n-1, l-1) - r1(1, -1)*r_prev(n-1, 1-l) - &
3205 & r1(-1, 1)*r_prev(1-n, l-1) + r1(-1, -1)*r_prev(1-n, 1-l)
3206 w = r1(1, 1)*r_prev(n+1, l-1) - r1(1, -1)*r_prev(n+1, 1-l) + &
3207 & r1(-1, 1)*r_prev(-n-1, l-1) - r1(-1, -1)*r_prev(-n-1, 1-l)
3208 r(n, l) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
3209 r(n, l) = r(n, l) / scal_uvw_m(l)
3211 u = r1(0, 1)*r_prev(-n, l-1) - r1(0, -1)*r_prev(-n, 1-l)
3212 v = r1(1, 1)*r_prev(1-n, l-1) - r1(1, -1)*r_prev(1-n, 1-l) + &
3213 & r1(-1, 1)*r_prev(n-1, l-1) - r1(-1, -1)*r_prev(n-1, 1-l)
3214 w = r1(1, 1)*r_prev(-n-1, l-1) - r1(1, -1)*r_prev(-n-1, 1-l) - &
3215 & r1(-1, 1)*r_prev(n+1, l-1) + r1(-1, -1)*r_prev(n+1, 1-l)
3216 r(-n, l) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
3217 r(-n, l) = r(-n, l) / scal_uvw_m(l)
3220 dst(ind+l) = alpha * dot_product(src(ind-l:ind+l), r(-l:l, l))
3222 v = r1(1, 1)*r_prev(l-1, 1-l) + r1(1, -1)*r_prev(l-1, l-1) - &
3223 & r1(-1, 1)*r_prev(1-l, 1-l) - r1(-1, -1)*r_prev(1-l, l-1)
3224 r(l, -l) = v * scal_v_n(l) / scal_uvw_m(l)
3227 v = r1(1, 1)*r_prev(1-l, 1-l) + r1(1, -1)*r_prev(1-l, l-1) + &
3228 & r1(-1, 1)*r_prev(l-1, 1-l) + r1(-1, -1)*r_prev(l-1, l-1)
3229 r(-l, -l) = v * scal_v_n(l) / scal_uvw_m(l)
3232 u = r1(0, 1)*r_prev(l-1, 1-l) + r1(0, -1)*r_prev(l-1, l-1)
3233 v = r1(1, 1)*r_prev(l-2, 1-l) + r1(1, -1)*r_prev(l-2, l-1) - &
3234 & r1(-1, 1)*r_prev(2-l, 1-l) - r1(-1, -1)*r_prev(2-l, l-1)
3235 r(l-1, -l) = (u*scal_u_n(l-1)+v*scal_v_n(l-1)) / scal_uvw_m(l)
3238 u = r1(0, 1)*r_prev(1-l, 1-l) + r1(0, -1)*r_prev(1-l, l-1)
3239 v = r1(1, 1)*r_prev(2-l, 1-l) + r1(1, -1)*r_prev(2-l, l-1) + &
3240 & r1(-1, 1)*r_prev(l-2, 1-l) + r1(-1, -1)*r_prev(l-2, l-1)
3241 r(1-l, -l) = (u*scal_u_n(l-1)+v*scal_v_n(l-1)) / scal_uvw_m(l)
3244 u = r1(0, 1)*r_prev(1, 1-l) + r1(0, -1)*r_prev(1, l-1)
3245 v = r1(1, 1)*r_prev(0, 1-l) + r1(1, -1)*r_prev(0, l-1)
3246 w = r1(1, 1)*r_prev(2, 1-l) + r1(1, -1)*r_prev(2, l-1) + &
3247 & r1(-1, 1)*r_prev(-2, 1-l) + r1(-1, -1)*r_prev(-2, l-1)
3248 r(1, -l) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
3249 r(1, -l) = r(1, -l) / scal_uvw_m(l)
3252 u = r1(0, 1)*r_prev(-1, 1-l) + r1(0, -1)*r_prev(-1, l-1)
3253 v = r1(-1, 1)*r_prev(0, 1-l) + r1(-1, -1)*r_prev(0, l-1)
3254 w = r1(1, 1)*r_prev(-2, 1-l) + r1(1, -1)*r_prev(-2, l-1) - &
3255 & r1(-1, 1)*r_prev(2, 1-l) - r1(-1, -1)*r_prev(2, l-1)
3256 r(-1, -l) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
3257 r(-1, -l) = r(-1, -l) / scal_uvw_m(l)
3260 u = r1(0, 1)*r_prev(0, 1-l) + r1(0, -1)*r_prev(0, l-1)
3261 v = r1(1, 1)*r_prev(1, 1-l) + r1(1, -1)*r_prev(1, l-1) + &
3262 & r1(-1, 1)*r_prev(-1, 1-l) + r1(-1, -1)*r_prev(-1, l-1)
3263 r(0, -l) = (u*scal_u_n(0) + v*scal_v_n(0)) / scal_uvw_m(l)
3267 u = r1(0, 1)*r_prev(n, 1-l) + r1(0, -1)*r_prev(n, l-1)
3268 v = r1(1, 1)*r_prev(n-1, 1-l) + r1(1, -1)*r_prev(n-1, l-1) - &
3269 & r1(-1, 1)*r_prev(1-n, 1-l) - r1(-1, -1)*r_prev(1-n, l-1)
3270 w = r1(1, 1)*r_prev(n+1, 1-l) + r1(1, -1)*r_prev(n+1, l-1) + &
3271 & r1(-1, 1)*r_prev(-n-1, 1-l) + r1(-1, -1)*r_prev(-n-1, l-1)
3272 r(n, -l) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
3273 r(n, -l) = r(n, -l) / scal_uvw_m(l)
3275 u = r1(0, 1)*r_prev(-n, 1-l) + r1(0, -1)*r_prev(-n, l-1)
3276 v = r1(1, 1)*r_prev(1-n, 1-l) + r1(1, -1)*r_prev(1-n, l-1) + &
3277 & r1(-1, 1)*r_prev(n-1, 1-l) + r1(-1, -1)*r_prev(n-1, l-1)
3278 w = r1(1, 1)*r_prev(-n-1, 1-l) + r1(1, -1)*r_prev(-n-1, l-1) - &
3279 & r1(-1, 1)*r_prev(n+1, 1-l) - r1(-1, -1)*r_prev(n+1, l-1)
3280 r(-n, -l) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
3281 r(-n, -l) = r(-n, -l) / scal_uvw_m(l)
3284 dst(ind-l) = alpha * dot_product(src(ind-l:ind+l), r(-l:l, -l))
3288 v = r1(1, 0)*r_prev(l-1, m) - r1(-1, 0)*r_prev(1-l, m)
3289 r(l, m) = v * scal_v_n(l) / scal_uvw_m(m)
3292 v = r1(1, 0)*r_prev(1-l, m) + r1(-1, 0)*r_prev(l-1, m)
3293 r(-l, m) = v * scal_v_n(l) / scal_uvw_m(m)
3296 u = r1(0, 0) * r_prev(l-1, m)
3297 v = r1(1, 0)*r_prev(l-2, m) - r1(-1, 0)*r_prev(2-l, m)
3298 r(l-1, m) = u*scal_u_n(l-1) + v*scal_v_n(l-1)
3299 r(l-1, m) = r(l-1, m) / scal_uvw_m(m)
3302 u = r1(0, 0) * r_prev(1-l, m)
3303 v = r1(1, 0)*r_prev(2-l, m) + r1(-1, 0)*r_prev(l-2, m)
3304 r(1-l, m) = u*scal_u_n(l-1) + v*scal_v_n(l-1)
3305 r(1-l, m) = r(1-l, m) / scal_uvw_m(m)
3308 u = r1(0, 0) * r_prev(0, m)
3309 v = r1(1, 0)*r_prev(1, m) + r1(-1, 0)*r_prev(-1, m)
3310 r(0, m) = u*scal_u_n(0) + v*scal_v_n(0)
3311 r(0, m) = r(0, m) / scal_uvw_m(m)
3314 u = r1(0, 0) * r_prev(1, m)
3315 v = r1(1, 0) * r_prev(0, m)
3316 w = r1(1, 0)*r_prev(2, m) + r1(-1, 0)*r_prev(-2, m)
3317 r(1, m) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
3318 r(1, m) = r(1, m) / scal_uvw_m(m)
3321 u = r1(0, 0) * r_prev(-1, m)
3322 v = r1(-1, 0) * r_prev(0, m)
3323 w = r1(1, 0)*r_prev(-2, m) - r1(-1, 0)*r_prev(2, m)
3324 r(-1, m) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
3325 r(-1, m) = r(-1, m) / scal_uvw_m(m)
3329 u = r1(0, 0) * r_prev(n, m)
3330 v = r1(1, 0)*r_prev(n-1, m) - r1(-1, 0)*r_prev(1-n, m)
3331 w = r1(1, 0)*r_prev(n+1, m) + r1(-1, 0)*r_prev(-1-n, m)
3332 r(n, m) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
3333 r(n, m) = r(n, m) / scal_uvw_m(m)
3335 u = r1(0, 0) * r_prev(-n, m)
3336 v = r1(1, 0)*r_prev(1-n, m) + r1(-1, 0)*r_prev(n-1, m)
3337 w = r1(1, 0)*r_prev(-n-1, m) - r1(-1, 0)*r_prev(n+1, m)
3338 r(-n, m) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
3339 r(-n, m) = r(-n, m) / scal_uvw_m(m)
3342 dst(ind+m) = alpha * dot_product(src(ind-l:ind+l), r(-l:l, m))
3348 dst(1) = beta*dst(1) + alpha*src(1)
3353 dst(2) = beta*dst(2) + alpha*dot_product(r1(-1:1,-1), src(2:4))
3354 dst(3) = beta*dst(3) + alpha*dot_product(r1(-1:1,0), src(2:4))
3355 dst(4) = beta*dst(4) + alpha*dot_product(r1(-1:1,1), src(2:4))
3362 r_prev(-p:p, -p:p) => work(1:l)
3364 r(-p:p, -p:p) => work(l+1:m)
3366 scal_uvw_m(-p:p) => work(m+1:l)
3368 scal_u_n(-p:p) => work(l+1:m)
3370 scal_v_n(-p:p) => work(m+1:l)
3372 scal_w_n(-p:p) => work(l+1:m)
3374 r(2, 2) = (r1(1, 1)*r1(1, 1) - r1(1, -1)*r1(1, -1) - &
3375 & r1(-1, 1)*r1(-1, 1) + r1(-1, -1)*r1(-1, -1)) / two
3376 r(2, 1) = r1(1, 1)*r1(1, 0) - r1(-1, 1)*r1(-1, 0)
3377 r(2, 0) = sqrt3 / two * (r1(1, 0)*r1(1, 0) - r1(-1, 0)*r1(-1, 0))
3378 r(2, -1) = r1(1, -1)*r1(1, 0) - r1(-1, -1)*r1(-1, 0)
3379 r(2, -2) = r1(1, 1)*r1(1, -1) - r1(-1, 1)*r1(-1, -1)
3380 r(1, 2) = r1(1, 1)*r1(0, 1) - r1(1, -1)*r1(0, -1)
3381 r(1, 1) = r1(1, 1)*r1(0, 0) + r1(1, 0)*r1(0, 1)
3382 r(1, 0) = sqrt3 * r1(1, 0) * r1(0, 0)
3383 r(1, -1) = r1(1, -1)*r1(0, 0) + r1(1, 0)*r1(0, -1)
3384 r(1, -2) = r1(1, 1)*r1(0, -1) + r1(1, -1)*r1(0, 1)
3385 r(0, 2) = sqrt3 / two * (r1(0, 1)*r1(0, 1) - r1(0, -1)*r1(0, -1))
3386 r(0, 1) = sqrt3 * r1(0, 1) * r1(0, 0)
3387 r(0, 0) = (three*r1(0, 0)*r1(0, 0)-one) / two
3388 r(0, -1) = sqrt3 * r1(0, -1) * r1(0, 0)
3389 r(0, -2) = sqrt3 * r1(0, 1) * r1(0, -1)
3390 r(-1, 2) = r1(-1, 1)*r1(0, 1) - r1(-1, -1)*r1(0, -1)
3391 r(-1, 1) = r1(-1, 1)*r1(0, 0) + r1(-1, 0)*r1(0, 1)
3392 r(-1, 0) = sqrt3 * r1(-1, 0) * r1(0, 0)
3393 r(-1, -1) = r1(-1, -1)*r1(0, 0) + r1(0, -1)*r1(-1, 0)
3394 r(-1, -2) = r1(-1, 1)*r1(0, -1) + r1(-1, -1)*r1(0, 1)
3395 r(-2, 2) = r1(1, 1)*r1(-1, 1) - r1(1, -1)*r1(-1, -1)
3396 r(-2, 1) = r1(1, 1)*r1(-1, 0) + r1(1, 0)*r1(-1, 1)
3397 r(-2, 0) = sqrt3 * r1(1, 0) * r1(-1, 0)
3398 r(-2, -1) = r1(1, -1)*r1(-1, 0) + r1(1, 0)*r1(-1, -1)
3399 r(-2, -2) = r1(1, 1)*r1(-1, -1) + r1(1, -1)*r1(-1, 1)
3401 dst(7+m) = beta*dst(7+m) + alpha*dot_product(r(-2:2, m), src(5:9))
3410 scal_uvw_m(0) = dble(l)
3412 u = sqrt(dble(l*l-m*m))
3417 u = sqrt(dble(u*(u-one)))
3420 scal_u_n(0) = dble(l)
3421 scal_v_n(0) = -sqrt(dble(l*(l-1))) / sqrt2
3424 u = sqrt(dble(l*l-n*n))
3428 v = sqrt(v*(v-one)) / two
3432 w = -sqrt(w*(w-one)) / two
3436 u = sqrt(dble(2*l-1))
3441 v = sqrt(dble((2*l-1)*(2*l-2))) / two
3444 v = sqrt(dble(2*l*(2*l-1))) / two
3447 scal_w_n(l-1) = zero
3450 scal_w_n(1-l) = zero
3453 v = r1(1, 1)*r_prev(l-1, l-1) - r1(1, -1)*r_prev(l-1, 1-l) - &
3454 & r1(-1, 1)*r_prev(1-l, l-1) + r1(-1, -1)*r_prev(1-l, 1-l)
3455 r(l, l) = v * scal_v_n(l) / scal_uvw_m(l)
3458 v = r1(1, 1)*r_prev(1-l, l-1) - r1(1, -1)*r_prev(1-l, 1-l) + &
3459 & r1(-1, 1)*r_prev(l-1, l-1) - r1(-1, -1)*r_prev(l-1, 1-l)
3460 r(-l, l) = v * scal_v_n(l) / scal_uvw_m(l)
3463 u = r1(0, 1)*r_prev(l-1, l-1) - r1(0, -1)*r_prev(l-1, 1-l)
3464 v = r1(1, 1)*r_prev(l-2, l-1) - r1(1, -1)*r_prev(l-2, 1-l) - &
3465 & r1(-1, 1)*r_prev(2-l, l-1) + r1(-1, -1)*r_prev(2-l, 1-l)
3466 r(l-1, l) = (u*scal_u_n(l-1)+v*scal_v_n(l-1)) / scal_uvw_m(l)
3469 u = r1(0, 1)*r_prev(1-l, l-1) - r1(0, -1)*r_prev(1-l, 1-l)
3470 v = r1(1, 1)*r_prev(2-l, l-1) - r1(1, -1)*r_prev(2-l, 1-l) + &
3471 & r1(-1, 1)*r_prev(l-2, l-1) - r1(-1, -1)*r_prev(l-2, 1-l)
3472 r(1-l, l) = (u*scal_u_n(l-1)+v*scal_v_n(l-1)) / scal_uvw_m(l)
3475 u = r1(0, 1)*r_prev(1, l-1) - r1(0, -1)*r_prev(1, 1-l)
3476 v = r1(1, 1)*r_prev(0, l-1) - r1(1, -1)*r_prev(0, 1-l)
3477 w = r1(1, 1)*r_prev(2, l-1) - r1(1, -1)*r_prev(2, 1-l) + &
3478 & r1(-1, 1)*r_prev(-2, l-1) - r1(-1, -1)*r_prev(-2, 1-l)
3479 r(1, l) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
3480 r(1, l) = r(1, l) / scal_uvw_m(l)
3483 u = r1(0, 1)*r_prev(-1, l-1) - r1(0, -1)*r_prev(-1, 1-l)
3484 v = r1(-1, 1)*r_prev(0, l-1) - r1(-1, -1)*r_prev(0, 1-l)
3485 w = r1(1, 1)*r_prev(-2, l-1) - r1(1, -1)*r_prev(-2, 1-l) - &
3486 & r1(-1, 1)*r_prev(2, l-1) + r1(-1, -1)*r_prev(2, 1-l)
3487 r(-1, l) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
3488 r(-1, l) = r(-1, l) / scal_uvw_m(l)
3491 u = r1(0, 1)*r_prev(0, l-1) - r1(0, -1)*r_prev(0, 1-l)
3492 v = r1(1, 1)*r_prev(1, l-1) - r1(1, -1)*r_prev(1, 1-l) + &
3493 & r1(-1, 1)*r_prev(-1, l-1) - r1(-1, -1)*r_prev(-1, 1-l)
3494 r(0, l) = (u*scal_u_n(0) + v*scal_v_n(0)) / scal_uvw_m(l)
3498 u = r1(0, 1)*r_prev(n, l-1) - r1(0, -1)*r_prev(n, 1-l)
3499 v = r1(1, 1)*r_prev(n-1, l-1) - r1(1, -1)*r_prev(n-1, 1-l) - &
3500 & r1(-1, 1)*r_prev(1-n, l-1) + r1(-1, -1)*r_prev(1-n, 1-l)
3501 w = r1(1, 1)*r_prev(n+1, l-1) - r1(1, -1)*r_prev(n+1, 1-l) + &
3502 & r1(-1, 1)*r_prev(-n-1, l-1) - r1(-1, -1)*r_prev(-n-1, 1-l)
3503 r(n, l) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
3504 r(n, l) = r(n, l) / scal_uvw_m(l)
3506 u = r1(0, 1)*r_prev(-n, l-1) - r1(0, -1)*r_prev(-n, 1-l)
3507 v = r1(1, 1)*r_prev(1-n, l-1) - r1(1, -1)*r_prev(1-n, 1-l) + &
3508 & r1(-1, 1)*r_prev(n-1, l-1) - r1(-1, -1)*r_prev(n-1, 1-l)
3509 w = r1(1, 1)*r_prev(-n-1, l-1) - r1(1, -1)*r_prev(-n-1, 1-l) - &
3510 & r1(-1, 1)*r_prev(n+1, l-1) + r1(-1, -1)*r_prev(n+1, 1-l)
3511 r(-n, l) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
3512 r(-n, l) = r(-n, l) / scal_uvw_m(l)
3515 dst(ind+l) = beta*dst(ind+l) + &
3516 & alpha*dot_product(src(ind-l:ind+l), r(-l:l, l))
3518 v = r1(1, 1)*r_prev(l-1, 1-l) + r1(1, -1)*r_prev(l-1, l-1) - &
3519 & r1(-1, 1)*r_prev(1-l, 1-l) - r1(-1, -1)*r_prev(1-l, l-1)
3520 r(l, -l) = v * scal_v_n(l) / scal_uvw_m(l)
3523 v = r1(1, 1)*r_prev(1-l, 1-l) + r1(1, -1)*r_prev(1-l, l-1) + &
3524 & r1(-1, 1)*r_prev(l-1, 1-l) + r1(-1, -1)*r_prev(l-1, l-1)
3525 r(-l, -l) = v * scal_v_n(l) / scal_uvw_m(l)
3528 u = r1(0, 1)*r_prev(l-1, 1-l) + r1(0, -1)*r_prev(l-1, l-1)
3529 v = r1(1, 1)*r_prev(l-2, 1-l) + r1(1, -1)*r_prev(l-2, l-1) - &
3530 & r1(-1, 1)*r_prev(2-l, 1-l) - r1(-1, -1)*r_prev(2-l, l-1)
3531 r(l-1, -l) = (u*scal_u_n(l-1)+v*scal_v_n(l-1)) / scal_uvw_m(l)
3534 u = r1(0, 1)*r_prev(1-l, 1-l) + r1(0, -1)*r_prev(1-l, l-1)
3535 v = r1(1, 1)*r_prev(2-l, 1-l) + r1(1, -1)*r_prev(2-l, l-1) + &
3536 & r1(-1, 1)*r_prev(l-2, 1-l) + r1(-1, -1)*r_prev(l-2, l-1)
3537 r(1-l, -l) = (u*scal_u_n(l-1)+v*scal_v_n(l-1)) / scal_uvw_m(l)
3540 u = r1(0, 1)*r_prev(1, 1-l) + r1(0, -1)*r_prev(1, l-1)
3541 v = r1(1, 1)*r_prev(0, 1-l) + r1(1, -1)*r_prev(0, l-1)
3542 w = r1(1, 1)*r_prev(2, 1-l) + r1(1, -1)*r_prev(2, l-1) + &
3543 & r1(-1, 1)*r_prev(-2, 1-l) + r1(-1, -1)*r_prev(-2, l-1)
3544 r(1, -l) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
3545 r(1, -l) = r(1, -l) / scal_uvw_m(l)
3548 u = r1(0, 1)*r_prev(-1, 1-l) + r1(0, -1)*r_prev(-1, l-1)
3549 v = r1(-1, 1)*r_prev(0, 1-l) + r1(-1, -1)*r_prev(0, l-1)
3550 w = r1(1, 1)*r_prev(-2, 1-l) + r1(1, -1)*r_prev(-2, l-1) - &
3551 & r1(-1, 1)*r_prev(2, 1-l) - r1(-1, -1)*r_prev(2, l-1)
3552 r(-1, -l) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
3553 r(-1, -l) = r(-1, -l) / scal_uvw_m(l)
3556 u = r1(0, 1)*r_prev(0, 1-l) + r1(0, -1)*r_prev(0, l-1)
3557 v = r1(1, 1)*r_prev(1, 1-l) + r1(1, -1)*r_prev(1, l-1) + &
3558 & r1(-1, 1)*r_prev(-1, 1-l) + r1(-1, -1)*r_prev(-1, l-1)
3559 r(0, -l) = (u*scal_u_n(0) + v*scal_v_n(0)) / scal_uvw_m(l)
3563 u = r1(0, 1)*r_prev(n, 1-l) + r1(0, -1)*r_prev(n, l-1)
3564 v = r1(1, 1)*r_prev(n-1, 1-l) + r1(1, -1)*r_prev(n-1, l-1) - &
3565 & r1(-1, 1)*r_prev(1-n, 1-l) - r1(-1, -1)*r_prev(1-n, l-1)
3566 w = r1(1, 1)*r_prev(n+1, 1-l) + r1(1, -1)*r_prev(n+1, l-1) + &
3567 & r1(-1, 1)*r_prev(-n-1, 1-l) + r1(-1, -1)*r_prev(-n-1, l-1)
3568 r(n, -l) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
3569 r(n, -l) = r(n, -l) / scal_uvw_m(l)
3571 u = r1(0, 1)*r_prev(-n, 1-l) + r1(0, -1)*r_prev(-n, l-1)
3572 v = r1(1, 1)*r_prev(1-n, 1-l) + r1(1, -1)*r_prev(1-n, l-1) + &
3573 & r1(-1, 1)*r_prev(n-1, 1-l) + r1(-1, -1)*r_prev(n-1, l-1)
3574 w = r1(1, 1)*r_prev(-n-1, 1-l) + r1(1, -1)*r_prev(-n-1, l-1) - &
3575 & r1(-1, 1)*r_prev(n+1, 1-l) - r1(-1, -1)*r_prev(n+1, l-1)
3576 r(-n, -l) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
3577 r(-n, -l) = r(-n, -l) / scal_uvw_m(l)
3580 dst(ind-l) = beta*dst(ind-l) + &
3581 & alpha*dot_product(src(ind-l:ind+l), r(-l:l, -l))
3585 v = r1(1, 0)*r_prev(l-1, m) - r1(-1, 0)*r_prev(1-l, m)
3586 r(l, m) = v * scal_v_n(l) / scal_uvw_m(m)
3589 v = r1(1, 0)*r_prev(1-l, m) + r1(-1, 0)*r_prev(l-1, m)
3590 r(-l, m) = v * scal_v_n(l) / scal_uvw_m(m)
3593 u = r1(0, 0) * r_prev(l-1, m)
3594 v = r1(1, 0)*r_prev(l-2, m) - r1(-1, 0)*r_prev(2-l, m)
3595 r(l-1, m) = u*scal_u_n(l-1) + v*scal_v_n(l-1)
3596 r(l-1, m) = r(l-1, m) / scal_uvw_m(m)
3599 u = r1(0, 0) * r_prev(1-l, m)
3600 v = r1(1, 0)*r_prev(2-l, m) + r1(-1, 0)*r_prev(l-2, m)
3601 r(1-l, m) = u*scal_u_n(l-1) + v*scal_v_n(l-1)
3602 r(1-l, m) = r(1-l, m) / scal_uvw_m(m)
3605 u = r1(0, 0) * r_prev(0, m)
3606 v = r1(1, 0)*r_prev(1, m) + r1(-1, 0)*r_prev(-1, m)
3607 r(0, m) = u*scal_u_n(0) + v*scal_v_n(0)
3608 r(0, m) = r(0, m) / scal_uvw_m(m)
3611 u = r1(0, 0) * r_prev(1, m)
3612 v = r1(1, 0) * r_prev(0, m)
3613 w = r1(1, 0)*r_prev(2, m) + r1(-1, 0)*r_prev(-2, m)
3614 r(1, m) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
3615 r(1, m) = r(1, m) / scal_uvw_m(m)
3618 u = r1(0, 0) * r_prev(-1, m)
3619 v = r1(-1, 0) * r_prev(0, m)
3620 w = r1(1, 0)*r_prev(-2, m) - r1(-1, 0)*r_prev(2, m)
3621 r(-1, m) = u*scal_u_n(1) + sqrt2*v*scal_v_n(1) + w*scal_w_n(1)
3622 r(-1, m) = r(-1, m) / scal_uvw_m(m)
3626 u = r1(0, 0) * r_prev(n, m)
3627 v = r1(1, 0)*r_prev(n-1, m) - r1(-1, 0)*r_prev(1-n, m)
3628 w = r1(1, 0)*r_prev(n+1, m) + r1(-1, 0)*r_prev(-1-n, m)
3629 r(n, m) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
3630 r(n, m) = r(n, m) / scal_uvw_m(m)
3632 u = r1(0, 0) * r_prev(-n, m)
3633 v = r1(1, 0)*r_prev(1-n, m) + r1(-1, 0)*r_prev(n-1, m)
3634 w = r1(1, 0)*r_prev(-n-1, m) - r1(-1, 0)*r_prev(n+1, m)
3635 r(-n, m) = u*scal_u_n(n) + v*scal_v_n(n) + w*scal_w_n(n)
3636 r(-n, m) = r(-n, m) / scal_uvw_m(m)
3639 dst(ind+m) = beta*dst(ind+m) + &
3640 & alpha*dot_product(src(ind-l:ind+l), r(-l:l, m))
3669 integer,
intent(in) :: p
3670 real(dp),
intent(in) :: vcos(p+1), vsin(p+1), alpha, src((p+1)*(p+1)), beta
3672 real(dp),
intent(inout) :: dst((p+1)*(p+1))
3700 integer,
intent(in) :: p
3701 real(dp),
intent(in) :: vcos(p+1), vsin(p+1), alpha, src((p+1)*(p+1)), beta
3703 real(dp),
intent(inout) :: dst((p+1)*(p+1))
3705 integer :: l, m, ind
3706 real(dp) :: v1, v2, v3, v4
3708 if (alpha .eq. zero)
then
3710 if (beta .eq. zero)
then
3720 if (beta .eq. zero)
then
3722 dst(1) = alpha*src(1)
3728 dst(ind) = alpha*src(ind)
3737 dst(ind+m) = alpha * (v1*v3-v2*v4)
3739 dst(ind-m) = alpha * (v1*v4+v2*v3)
3744 dst(1) = beta*dst(1) + alpha*src(1)
3750 dst(ind) = beta*dst(ind) + alpha*src(ind)
3759 dst(ind+m) = beta*dst(ind+m) + alpha*(v1*v3-v2*v4)
3761 dst(ind-m) = beta*dst(ind-m) + alpha*(v1*v4+v2*v3)
3790 integer,
intent(in) :: p
3791 real(dp),
intent(in) :: vcos(p+1), vsin(p+1), alpha, src((p+1)*(p+1)), beta
3793 real(dp),
intent(inout) :: dst((p+1)*(p+1))
3821 integer,
intent(in) :: p
3822 real(dp),
intent(in) :: vcos(p+1), vsin(p+1), alpha, src((p+1)*(p+1)), beta
3824 real(dp),
intent(inout) :: dst((p+1)*(p+1))
3826 integer :: l, m, ind
3827 real(dp) :: v1, v2, v3, v4
3829 if (alpha .eq. zero)
then
3831 if (beta .eq. zero)
then
3841 if (beta .eq. zero)
then
3843 dst(1) = alpha*src(1)
3848 dst(ind) = alpha*src(ind)
3856 dst(ind+m) = alpha * (v1*v3+v2*v4)
3858 dst(ind-m) = alpha * (v2*v3-v1*v4)
3863 dst(1) = beta*dst(1) + alpha*src(1)
3868 dst(ind) = beta*dst(ind) + alpha*src(ind)
3876 dst(ind+m) = beta*dst(ind+m) + alpha*(v1*v3+v2*v4)
3878 dst(ind-m) = beta*dst(ind-m) + alpha*(v2*v3-v1*v4)
3918 integer,
intent(in) :: p
3919 real(dp),
intent(in) :: ctheta, stheta, alpha, src((p+1)**2), beta
3921 real(dp),
intent(inout) :: dst((p+1)**2)
3923 real(dp) :: work(2*(2*p+1)*(2*p+3)+p)
3944 integer,
intent(in) :: p
3945 real(dp),
intent(in) :: ctheta, stheta, alpha, src((p+1)**2), beta
3947 real(dp),
intent(out) :: dst((p+1)*(p+1))
3949 real(dp),
intent(out),
target :: work(4*p*p+13*p+4)
3951 real(dp) :: u, v, w, fl, fl2, tmp1, tmp2, vu(2), vv(2), vw(2), &
3952 ctheta2, stheta2, cstheta
3953 integer :: l, m, n, ind
3955 real(dp),
pointer :: r(:, :, :), r_prev(:, :, :), scal_uvw_m(:), &
3956 & scal_u_n(:), scal_v_n(:), scal_w_n(:), r_swap(:, :, :), vsqr(:)
3963 if (alpha .eq. zero)
then
3965 if (beta .eq. zero)
then
3975 if (beta .eq. zero)
then
3978 dst(1) = alpha * src(1)
3983 dst(2) = alpha * src(2)
3984 dst(3) = alpha * (src(3)*ctheta - src(4)*stheta)
3985 dst(4) = alpha * (src(3)*stheta + src(4)*ctheta)
3990 l = 2 * (p+1) * (p+1)
3991 r(1:2, 0:p, 0:p) => work(1:l)
3993 r_prev(1:2, 0:p, 0:p) => work(l+1:m)
3995 scal_uvw_m(0:p) => work(m+1:l)
3997 scal_u_n(0:p-1) => work(l+1:m)
3999 scal_v_n(0:p) => work(m+1:l)
4001 scal_w_n(1:p-2) => work(l+1:m)
4003 vsqr(1:p) => work(m+1:l)
4005 ctheta2 = ctheta * ctheta
4006 cstheta = ctheta * stheta
4007 stheta2 = stheta * stheta
4008 r(1, 2, 2) = (ctheta2 + one) / two
4009 r(1, 1, 2) = cstheta
4010 r(1, 0, 2) = sqrt3 / two * stheta2
4011 dst(9) = alpha * (src(9)*r(1, 2, 2) + src(8)*r(1, 1, 2) + &
4012 & src(7)*r(1, 0, 2))
4013 r(1, 2, 1) = -cstheta
4014 r(1, 1, 1) = ctheta2 - stheta2
4015 r(1, 0, 1) = sqrt3 * cstheta
4016 dst(8) = alpha * (src(9)*r(1, 2, 1) + src(8)*r(1, 1, 1) + &
4017 & src(7)*r(1, 0, 1))
4018 r(1, 2, 0) = sqrt3 / two * stheta2
4019 r(1, 1, 0) = -sqrt3 * cstheta
4020 r(1, 0, 0) = (three*ctheta2-one) / two
4021 dst(7) = alpha * (src(9)*r(1, 2, 0) + src(8)*r(1, 1, 0) + &
4022 & src(7)*r(1, 0, 0))
4025 r(2, 2, 1) = -stheta
4026 dst(6) = alpha * (src(6)*r(2, 1, 1) + src(5)*r(2, 2, 1))
4029 dst(5) = alpha * (src(6)*r(2, 1, 2) + src(5)*r(2, 2, 2))
4042 scal_uvw_m(0) = one / fl
4045 u = sqrt(fl2 - vsqr(m))
4046 scal_uvw_m(m) = one / u
4049 u = sqrt(dble(u*(u-one)))
4050 scal_uvw_m(l) = one / u
4051 scal_u_n(0) = dble(l)
4052 scal_v_n(0) = -sqrt(dble(l*(l-1))) / sqrt2
4055 u = sqrt(fl2-vsqr(n))
4061 v = sqrt(v*v-v) / two
4064 w = -sqrt(w*w-w) / two
4067 u = sqrt(dble(2*l-1))
4069 v = sqrt(dble((2*l-1)*(2*l-2))) / two
4071 v = sqrt(dble(2*l*(2*l-1))) / two
4075 vv = ctheta*r_prev(:, l-1, l-1) + r_prev(2:1:-1, l-1, l-1)
4076 r(:, l, l) = vv * scal_v_n(l) * scal_uvw_m(l)
4077 tmp1 = src(ind+l) * r(1, l, l)
4078 tmp2 = src(ind-l) * r(2, l, l)
4080 vu = stheta * r_prev(:, l-1, l-1)
4081 vv = ctheta*r_prev(:, l-2, l-1) + r_prev(2:1:-1, l-2, l-1)
4082 r(:, l-1, l) = (vu*scal_u_n(l-1)+vv*scal_v_n(l-1)) * scal_uvw_m(l)
4083 tmp1 = tmp1 + src(ind+l-1)*r(1, l-1, l)
4084 tmp2 = tmp2 + src(ind-l+1)*r(2, l-1, l)
4086 vu = stheta * r_prev(:, 1, l-1)
4087 vv(1) = ctheta * r_prev(1, 0, l-1)
4088 vv(2) = r_prev(1, 0, l-1)
4089 vw = ctheta*r_prev(:, 2, l-1) - r_prev(2:1:-1, 2, l-1)
4090 r(:, 1, l) = vu*scal_u_n(1) + vw*scal_w_n(1) + sqrt2*scal_v_n(1)*vv
4091 r(:, 1, l) = r(:, 1, l) * scal_uvw_m(l)
4092 tmp1 = tmp1 + src(ind+1)*r(1, 1, l)
4093 tmp2 = tmp2 + src(ind-1)*r(2, 1, l)
4095 u = stheta * r_prev(1, 0, l-1)
4096 v = ctheta*r_prev(1, 1, l-1) - r_prev(2, 1, l-1)
4097 r(1, 0, l) = (u*scal_u_n(0) + v*scal_v_n(0)) * scal_uvw_m(l)
4098 tmp1 = tmp1 + src(ind)*r(1, 0, l)
4102 vu = stheta * r_prev(:, n, l-1)
4103 vv = ctheta*r_prev(:, n-1, l-1) + r_prev(2:1:-1, n-1, l-1)
4104 vw = ctheta*r_prev(:, n+1, l-1) - r_prev(2:1:-1, n+1, l-1)
4105 vu = vu*scal_u_n(n) + vv*scal_v_n(n) + vw*scal_w_n(n)
4106 r(:, n, l) = vu * scal_uvw_m(l)
4107 tmp1 = tmp1 + src(ind+n)*r(1, n, l)
4108 tmp2 = tmp2 + src(ind-n)*r(2, n, l)
4110 dst(ind+l) = alpha * tmp1
4111 dst(ind-l) = alpha * tmp2
4114 v = -stheta * r_prev(1, l-1, 0)
4115 u = scal_v_n(l) * scal_uvw_m(0)
4117 tmp1 = src(ind+l) * r(1, l, 0)
4119 u = ctheta * r_prev(1, l-1, 0)
4120 v = -stheta * r_prev(1, l-2, 0)
4121 w = u*scal_u_n(l-1) + v*scal_v_n(l-1)
4122 r(1, l-1, 0) = w * scal_uvw_m(0)
4123 tmp1 = tmp1 + src(ind+l-1)*r(1, l-1, 0)
4125 u = ctheta * r_prev(1, 0, 0)
4126 v = -stheta * r_prev(1, 1, 0)
4127 w = u*scal_u_n(0) + v*scal_v_n(0)
4128 r(1, 0, 0) = w * scal_uvw_m(0)
4129 tmp1 = tmp1 + src(ind)*r(1, 0, 0)
4131 v = sqrt2*scal_v_n(1)*r_prev(1, 0, 0) + &
4132 & scal_w_n(1)*r_prev(1, 2, 0)
4133 u = ctheta * r_prev(1, 1, 0)
4134 w = scal_u_n(1)*u - stheta*v
4135 r(1, 1, 0) = w * scal_uvw_m(0)
4136 tmp1 = tmp1 + src(ind+1)*r(1, 1, 0)
4140 v = scal_v_n(n)*r_prev(1, n-1, 0) + &
4141 & scal_w_n(n)*r_prev(1, n+1, 0)
4142 u = ctheta * r_prev(1, n, 0)
4143 w = scal_u_n(n)*u - stheta*v
4144 r(1, n, 0) = w * scal_uvw_m(0)
4145 tmp1 = tmp1 + src(ind+n)*r(1, n, 0)
4147 dst(ind) = alpha * tmp1
4152 vv = -stheta * r_prev(:, l-1, m)
4153 u = scal_v_n(l) * scal_uvw_m(m)
4155 tmp1 = src(ind+l) * r(1, l, m)
4156 tmp2 = src(ind-l) * r(2, l, m)
4158 vu = ctheta * r_prev(:, l-1, m)
4159 vv = -stheta * r_prev(:, l-2, m)
4160 vw = vu*scal_u_n(l-1) + vv*scal_v_n(l-1)
4161 r(:, l-1, m) = vw * scal_uvw_m(m)
4162 tmp1 = tmp1 + src(ind+l-1)*r(1, l-1, m)
4163 tmp2 = tmp2 + src(ind-l+1)*r(2, l-1, m)
4165 u = ctheta * r_prev(1, 0, m)
4166 v = -stheta * r_prev(1, 1, m)
4167 w = u*scal_u_n(0) + v*scal_v_n(0)
4168 r(1, 0, m) = w * scal_uvw_m(m)
4169 tmp1 = tmp1 + src(ind)*r(1, 0, m)
4171 v = sqrt2*scal_v_n(1)*r_prev(1, 0, m) + &
4172 & scal_w_n(1)*r_prev(1, 2, m)
4173 u = ctheta * r_prev(1, 1, m)
4174 w = scal_u_n(1)*u - stheta*v
4175 r(1, 1, m) = w * scal_uvw_m(m)
4176 tmp1 = tmp1 + src(ind+1)*r(1, 1, m)
4178 u = ctheta * r_prev(2, 1, m)
4179 w = -stheta * r_prev(2, 2, m)
4180 v = u*scal_u_n(1) + w*scal_w_n(1)
4181 r(2, 1, m) = v * scal_uvw_m(m)
4182 tmp2 = tmp2 + src(ind-1)*r(2, 1, m)
4186 vv = scal_v_n(n)*r_prev(:, n-1, m) + &
4187 & scal_w_n(n)*r_prev(:, n+1, m)
4188 vu = ctheta * r_prev(:, n, m)
4189 vw = scal_u_n(n)*vu - stheta*vv
4190 r(:, n, m) = vw * scal_uvw_m(m)
4191 tmp1 = tmp1 + src(ind+n)*r(1, n, m)
4192 tmp2 = tmp2 + src(ind-n)*r(2, n, m)
4194 dst(ind+m) = alpha * tmp1
4195 dst(ind-m) = alpha * tmp2
4199 stop
"Not Implemented"
4226 & src_m, beta, dst_m)
4228 integer,
intent(in) :: p
4229 real(dp),
intent(in) :: z, src_r, dst_r, vscales((p+1)*(p+1)), &
4230 & vcnk((2*p+1)*(p+1)), alpha, src_m((p+1)*(p+1)), beta
4232 real(dp),
intent(inout) :: dst_m((p+1)*(p+1))
4234 real(dp) :: work(2*(p+1))
4237 & src_m, beta, dst_m, work)
4264 & src_m, beta, dst_m, work)
4266 integer,
intent(in) :: p
4267 real(dp),
intent(in) :: z, src_r, dst_r, vscales((p+1)*(p+1)), &
4268 & vcnk((2*p+1)*(p+1)), alpha, src_m((p+1)*(p+1)), beta
4270 real(dp),
intent(inout) :: dst_m((p+1)*(p+1))
4272 real(dp),
intent(out),
target :: work(2*(p+1))
4274 real(dp) :: r1, r2, tmp1, tmp2, tmp3, res1, res2, pow_r1
4275 integer :: j, k, n, indj, indjn, indjk1, indjk2
4277 real(dp),
pointer :: pow_r2(:)
4279 if (alpha .eq. zero)
then
4280 if (beta .eq. zero)
then
4283 dst_m = beta * dst_m
4291 pow_r2(1:p+1) => work(1:p+1)
4300 pow_r2(j) = pow_r2(j-1) * r2
4304 if (beta .eq. zero)
then
4309 tmp1 = alpha * pow_r1 * vscales(indj)
4310 pow_r1 = pow_r1 * r1
4314 indjk1 = j*(j+1)/2 + 1
4317 indjn = (j-n)**2 + (j-n) + 1
4318 tmp3 = pow_r2(n+1) / &
4319 & vscales(indjn) * vcnk(indjk1+n)**2
4320 res1 = res1 + tmp3*src_m(indjn)
4322 dst_m(indj) = tmp2 * res1
4329 indjk1 = (j-k)*(j-k+1)/2 + 1
4330 indjk2 = (j+k)*(j+k+1)/2 + 1
4333 indjn = (j-n)**2 + (j-n) + 1
4334 tmp3 = pow_r2(n+1) / &
4335 & vscales(indjn) * vcnk(indjk1+n) * &
4337 res1 = res1 + tmp3*src_m(indjn+k)
4338 res2 = res2 + tmp3*src_m(indjn-k)
4340 dst_m(indj+k) = tmp2 * res1
4341 dst_m(indj-k) = tmp2 * res2
4350 tmp1 = alpha * pow_r1 * vscales(indj)
4351 pow_r1 = pow_r1 * r1
4355 indjk1 = j * (j+1) /2 + 1
4358 indjn = (j-n)**2 + (j-n) + 1
4359 tmp3 = pow_r2(n+1) / &
4360 & vscales(indjn) * vcnk(indjk1+n)**2
4361 res1 = res1 + tmp3*src_m(indjn)
4363 dst_m(indj) = beta*dst_m(indj) + tmp2*res1
4370 indjk1 = (j-k)*(j-k+1)/2 + 1
4371 indjk2 = (j+k)*(j+k+1)/2 + 1
4374 indjn = (j-n)**2 + (j-n) + 1
4375 tmp3 = pow_r2(n+1) / &
4376 & vscales(indjn) * vcnk(indjk1+n) * &
4378 res1 = res1 + tmp3*src_m(indjn+k)
4379 res2 = res2 + tmp3*src_m(indjn-k)
4381 dst_m(indj+k) = beta*dst_m(indj+k) + tmp2*res1
4382 dst_m(indj-k) = beta*dst_m(indj-k) + tmp2*res2
4389 if (beta .eq. zero)
then
4394 do k = indj-j, indj+j
4395 dst_m(k) = tmp1 * src_m(k)
4405 do k = indj-j, indj+j
4406 dst_m(k) = beta*dst_m(k) + tmp1*src_m(k)
4436 & src_m, beta, dst_m)
4438 integer,
intent(in) :: p
4439 real(dp),
intent(in) :: z, src_sk(p+1), dst_sk(p+1), vscales((p+1)*(p+1)), &
4440 & alpha, src_m((p+1)*(p+1)), beta
4442 real(dp),
intent(inout) :: dst_m((p+1)*(p+1))
4444 real(dp) :: work(2*p+1)
4445 complex(dp) :: work_complex(2*p+1)
4448 & alpha, src_m, beta, dst_m, work, work_complex)
4475 & alpha, src_m, beta, dst_m, work, work_complex)
4478 integer,
intent(in) :: p
4479 real(dp),
intent(in) :: z, src_sk(p+1), dst_sk(p+1), vscales((p+1)*(p+1)), &
4480 & alpha, src_m((p+1)*(p+1)), beta
4482 real(dp),
intent(inout) :: dst_m((p+1)*(p+1))
4484 real(dp),
intent(out),
target :: work(2*p+1)
4485 complex(dp),
intent(out) :: work_complex(2*p+1)
4487 real(dp) :: r1, tmp1, tmp2, tmp3, res1, res2
4488 integer :: j, k, n, indj, indn, l, NZ, ierr
4490 real(dp) :: fact(2*p+1)
4491 complex(dp) :: complex_argument
4493 if (alpha .eq. zero)
then
4494 if (beta .eq. zero)
then
4497 dst_m = beta * dst_m
4504 fact(j+1) = dble(j) * fact(j)
4509 complex_argument = z
4510 call cbesi(complex_argument, pt5, 1, 1, work_complex(1), nz, ierr)
4511 work(1) = sqrt(pi/(2*z)) * real(work_complex(1))
4513 call cbesi(complex_argument, 1.5d0, 1, 2*p, &
4514 & work_complex(2:2*p+1), nz, ierr)
4515 work(2:2*p+1) = sqrt(pi/(2*z)) * real(work_complex(2:2*p+1))
4519 if (beta .eq. zero)
then
4530 tmp1 = vscales(indn) * &
4531 & dble(2*j+1) * fact(j+1) * fact(n+1) / &
4532 & src_sk(n+1) / vscales(indj) * dst_sk(j+1)
4535 tmp2 = (two**(-l)) / &
4536 & fact(l+1)*fact(2*l+1)/ &
4537 & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
4538 & (work(j+n-l+1) / (z**l))
4542 res1 = res1 + tmp3*src_m(indn)
4552 tmp1 = vscales(indn+k) * &
4553 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
4554 & src_sk(n+1) / vscales(indj+k) * dst_sk(j+1)
4557 tmp2 = (two**(-l)) / &
4558 & fact(l+k+1)*fact(2*l+1)/ &
4559 & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
4560 & (work(j+n-l+1) / (z**l))
4564 res1 = res1 + tmp3*src_m(indn+k)
4565 res2 = res2 + tmp3*src_m(indn-k)
4567 dst_m(indj+k) = res1
4568 dst_m(indj-k) = res2
4582 tmp1 = vscales(indn) * &
4583 & dble(2*j+1) * fact(j+1) * fact(n+1) / &
4584 & src_sk(n+1) / vscales(indj) * dst_sk(j+1)
4587 tmp2 = (two**(-l)) / &
4588 & fact(l+1)*fact(2*l+1)/ &
4589 & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
4590 & (work(j+n-l+1) / (z**l))
4594 res1 = res1 + tmp3*src_m(indn)
4596 dst_m(indj) = beta*dst_m(indj) + res1
4604 tmp1 = vscales(indn+k) * &
4605 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
4606 & src_sk(n+1) / vscales(indj+k) * dst_sk(j+1)
4609 tmp2 = (two**(-l)) / &
4610 & fact(l+k+1)*fact(2*l+1)/ &
4611 & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
4612 & (work(j+n-l+1) / (z**l))
4616 res1 = res1 + tmp3*src_m(indn+k)
4617 res2 = res2 + tmp3*src_m(indn-k)
4619 dst_m(indj+k) = beta*dst_m(indj+k) + res1
4620 dst_m(indj-k) = beta*dst_m(indj-k) + res2
4628 if (beta .eq. zero)
then
4633 do k = indj-j, indj+j
4634 dst_m(k) = tmp1 * src_m(k)
4644 do k = indj-j, indj+j
4645 dst_m(k) = beta*dst_m(k) + tmp1*src_m(k)
4675 & src_m, beta, dst_m)
4677 integer,
intent(in) :: p
4678 real(dp),
intent(in) :: z, src_sk(p+1), dst_sk(p+1), vscales((p+1)*(p+1)), &
4679 & alpha, src_m((p+1)*(p+1)), beta
4681 real(dp),
intent(inout) :: dst_m((p+1)*(p+1))
4683 real(dp) :: work(2*p+1)
4684 complex(dp) :: work_complex(2*p+1)
4687 & alpha, src_m, beta, dst_m, work, work_complex)
4714 & alpha, src_m, beta, dst_m, work, work_complex)
4717 integer,
intent(in) :: p
4718 real(dp),
intent(in) :: z, src_sk(p+1), dst_sk(p+1), vscales((p+1)*(p+1)), &
4719 & alpha, src_m((p+1)*(p+1)), beta
4721 real(dp),
intent(inout) :: dst_m((p+1)*(p+1))
4723 real(dp),
intent(out),
target :: work(2*p+1)
4724 complex(dp),
intent(out) :: work_complex(2*p+1)
4726 real(dp) :: r1, tmp1, tmp2, tmp3, res1, res2
4727 integer :: j, k, n, indj, indn, l, NZ, ierr
4729 real(dp) :: fact(2*p+1)
4730 complex(dp) :: complex_argument
4732 if (alpha .eq. zero)
then
4733 if (beta .eq. zero)
then
4736 dst_m = beta * dst_m
4743 fact(j+1) = dble(j) * fact(j)
4748 complex_argument = z
4749 call cbesi(complex_argument, pt5, 1, 1, work_complex(1), nz, ierr)
4750 work(1) = sqrt(pi/(2*z)) * real(work_complex(1))
4752 call cbesi(complex_argument, 1.5d0, 1, 2*p, &
4753 & work_complex(2:2*p+1), nz, ierr)
4754 work(2:2*p+1) = sqrt(pi/(2*z)) * real(work_complex(2:2*p+1))
4758 if (beta .eq. zero)
then
4769 tmp1 = vscales(indn) * &
4770 & dble(2*j+1) * fact(j+1) * fact(n+1) / &
4771 & src_sk(n+1) / vscales(indj) * dst_sk(j+1)
4774 tmp2 = (two**(-l)) / &
4775 & fact(l+1)*fact(2*l+1)/ &
4776 & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
4777 & (work(j+n-l+1) / (z**l))
4781 res1 = res1 + tmp3*src_m(indj)
4791 tmp1 = vscales(indn+k) * &
4792 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
4793 & src_sk(n+1) / vscales(indj+k) * dst_sk(j+1)
4796 tmp2 = (two**(-l)) / &
4797 & fact(l+k+1)*fact(2*l+1)/ &
4798 & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
4799 & (work(j+n-l+1) / (z**l))
4803 res1 = res1 + tmp3*src_m(indj+k)
4804 res2 = res2 + tmp3*src_m(indj-k)
4806 dst_m(indn+k) = res1
4807 dst_m(indn-k) = res2
4821 tmp1 = vscales(indn) * &
4822 & dble(2*j+1) * fact(j+1) * fact(n+1) / &
4823 & src_sk(n+1) / vscales(indj) * dst_sk(j+1)
4826 tmp2 = (two**(-l)) / &
4827 & fact(l+1)*fact(2*l+1)/ &
4828 & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
4829 & (work(j+n-l+1) / (z**l))
4833 res1 = res1 + tmp3*src_m(indj)
4835 dst_m(indn) = beta*dst_m(indn) + res1
4843 tmp1 = vscales(indn+k) * &
4844 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
4845 & src_sk(n+1) / vscales(indj+k) * dst_sk(j+1)
4848 tmp2 = (two**(-l)) / &
4849 & fact(l+k+1)*fact(2*l+1)/ &
4850 & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
4851 & (work(j+n-l+1) / (z**l))
4855 res1 = res1 + tmp3*src_m(indj+k)
4856 res2 = res2 + tmp3*src_m(indj-k)
4858 dst_m(indn+k) = beta*dst_m(indn+k) + res1
4859 dst_m(indn-k) = beta*dst_m(indn-k) + res2
4867 if (beta .eq. zero)
then
4872 do k = indj-j, indj+j
4873 dst_m(k) = tmp1 * src_m(k)
4883 do k = indj-j, indj+j
4884 dst_m(k) = beta*dst_m(k) + tmp1*src_m(k)
4916 & alpha, src_m, beta, dst_m)
4919 integer,
intent(in) :: p
4920 real(dp),
intent(in) :: src_sk(p+2), vscales((p+2)*(p+2)), &
4921 & alpha, src_m((p+1)*(p+1)), beta
4923 real(dp),
intent(inout) :: dst_m((p+2)*(p+2))
4926 real(dp) :: tmp1, tmp2, tmp3, res1, res2
4927 integer :: j, k, n, indj, indn, l
4929 real(dp) :: fact(2*p+1), fact2(p+2)
4931 if (alpha .eq. zero)
then
4932 if (beta .eq. zero)
then
4935 dst_m = beta * dst_m
4942 fact(j+1) = dble(j) * fact(j)
4946 fact2(j+1) = dble(2*j+1) * fact2(j)
4951 if (beta .eq. zero)
then
4960 dst_m(1) = vscales(3) / src_sk(2) / vscales(1) * src_sk(1) / fact2(2) * &
4973 tmp1 = vscales(indn) * &
4974 & dble(2*j+1) * fact(j+1) * fact(n+1) / &
4975 & src_sk(n+1) / vscales(indj) * src_sk(j+1) * alpha
4979 tmp2 = (two**(-l)) / &
4980 & fact(l+1)*fact(2*l+1)/ &
4981 & fact(l+1)/fact(l+1)/ &
4985 res1 = res1 + tmp3*src_m(indn)
4990 tmp1 = vscales(indn) * &
4991 & dble(2*j+1) * fact(j+1) * fact(n+1) / &
4992 & src_sk(n+1) / vscales(indj) * src_sk(j+1) * alpha
4996 tmp2 = (two**(-l)) / &
4997 & fact(l+1)*fact(2*l+1)/ &
4998 & fact(l+1)/fact(l+1)/ &
5002 res1 = res1 + tmp3*src_m(indn)
5012 tmp1 = vscales(indn+k) * &
5013 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
5014 & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
5017 tmp2 = (two**(-l)) / &
5018 & fact(l+k+1)*fact(2*l+1)/ &
5019 & fact(l+1)/fact(l-k+1)/ &
5023 res1 = res1 + tmp3*src_m(indn+k)
5024 res2 = res2 + tmp3*src_m(indn-k)
5029 tmp1 = vscales(indn+k) * &
5030 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
5031 & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
5034 tmp2 = (two**(-l)) / &
5035 & fact(l+k+1)*fact(2*l+1)/ &
5036 & fact(l+1)/fact(l-k+1)/ &
5040 res1 = res1 + tmp3*src_m(indn+k)
5041 res2 = res2 + tmp3*src_m(indn-k)
5042 dst_m(indj+k) = res1
5043 dst_m(indj-k) = res2
5053 tmp1 = vscales(indn+k) * &
5054 & dble(2*j+1) * fact(n+k+1) / &
5055 & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
5058 tmp2 = (two**(-l)) / &
5059 & fact(l+k+1)*fact(2*l+1)/ &
5064 res1 = res1 + tmp3*src_m(indn+k)
5065 res2 = res2 + tmp3*src_m(indn-k)
5066 dst_m(indj+k) = res1
5067 dst_m(indj-k) = res2
5079 tmp1 = vscales(indn) * &
5080 & dble(2*j+1) * fact(j+1) * fact(n+1) / &
5081 & src_sk(n+1) / vscales(indj) * src_sk(j+1) * alpha
5084 tmp2 = (two**(-l)) / &
5085 & fact(l+1)*fact(2*l+1)/ &
5086 & fact(l+1)/fact(l+1)/ &
5090 res1 = res1 + tmp3*src_m(indn)
5100 tmp1 = vscales(indn+k) * &
5101 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
5102 & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
5105 tmp2 = (two**(-l)) / &
5106 & fact(l+k+1)*fact(2*l+1)/ &
5107 & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)/ &
5111 res1 = res1 + tmp3*src_m(indn+k)
5112 res2 = res2 + tmp3*src_m(indn-k)
5113 dst_m(indj+k) = res1
5114 dst_m(indj-k) = res2
5127 tmp1 = vscales(indn) * &
5128 & dble(2*j+1) * fact(j+1) * fact(n+1) / &
5129 & src_sk(n+1) / vscales(indj) * src_sk(j+1) * alpha
5132 tmp2 = (two**(-l)) / &
5133 & fact(l+1)*fact(2*l+1)/ &
5134 & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)/ &
5138 res1 = res1 + tmp3*src_m(indn)
5148 tmp1 = vscales(indn+k) * &
5149 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
5150 & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
5153 tmp2 = (two**(-l)) / &
5154 & fact(l+k+1)*fact(2*l+1)/ &
5155 & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)/ &
5159 res1 = res1 + tmp3*src_m(indn+k)
5160 res2 = res2 + tmp3*src_m(indn-k)
5161 dst_m(indj+k) = res1
5162 dst_m(indj-k) = res2
5173 dst_m(1) = vscales(3) / src_sk(2) / vscales(1) * src_sk(1) / fact2(2) * &
5174 & src_m(3) * alpha + dst_m(1)*beta
5186 tmp1 = vscales(indn) * &
5187 & dble(2*j+1) * fact(j+1) * fact(n+1) / &
5188 & src_sk(n+1) / vscales(indj) * src_sk(j+1) * alpha
5192 tmp2 = (two**(-l)) / &
5193 & fact(l+1)*fact(2*l+1)/ &
5194 & fact(l+1)/fact(l+1)/ &
5198 res1 = res1 + tmp3*src_m(indn)
5203 tmp1 = vscales(indn) * &
5204 & dble(2*j+1) * fact(j+1) * fact(n+1) / &
5205 & src_sk(n+1) / vscales(indj) * src_sk(j+1) * alpha
5209 tmp2 = (two**(-l)) / &
5210 & fact(l+1)*fact(2*l+1)/ &
5211 & fact(l+1)/fact(l+1)/ &
5215 res1 = res1 + tmp3*src_m(indn)
5216 dst_m(indj) = beta*dst_m(indj) + res1
5225 tmp1 = vscales(indn+k) * &
5226 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
5227 & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
5230 tmp2 = (two**(-l)) / &
5231 & fact(l+k+1)*fact(2*l+1)/ &
5232 & fact(l+1)/fact(l-k+1)/ &
5236 res1 = res1 + tmp3*src_m(indn+k)
5237 res2 = res2 + tmp3*src_m(indn-k)
5242 tmp1 = vscales(indn+k) * &
5243 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
5244 & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
5247 tmp2 = (two**(-l)) / &
5248 & fact(l+k+1)*fact(2*l+1)/ &
5249 & fact(l+1)/fact(l-k+1)/ &
5253 res1 = res1 + tmp3*src_m(indn+k)
5254 res2 = res2 + tmp3*src_m(indn-k)
5255 dst_m(indj+k) = beta*dst_m(indj+k) + res1
5256 dst_m(indj-k) = beta*dst_m(indj-k) + res2
5266 tmp1 = vscales(indn+k) * &
5267 & dble(2*j+1) * fact(n+k+1) / &
5268 & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
5271 tmp2 = (two**(-l)) / &
5272 & fact(l+k+1)*fact(2*l+1)/ &
5277 res1 = res1 + tmp3*src_m(indn+k)
5278 res2 = res2 + tmp3*src_m(indn-k)
5279 dst_m(indj+k) = beta*dst_m(indj+k) + res1
5280 dst_m(indj-k) = beta*dst_m(indj-k) + res2
5292 tmp1 = vscales(indn) * &
5293 & dble(2*j+1) * fact(j+1) * fact(n+1) / &
5294 & src_sk(n+1) / vscales(indj) * src_sk(j+1) * alpha
5297 tmp2 = (two**(-l)) / &
5298 & fact(l+1)*fact(2*l+1)/ &
5299 & fact(l+1)/fact(l+1)/ &
5303 res1 = res1 + tmp3*src_m(indn)
5304 dst_m(indj) = beta*dst_m(indj) + res1
5313 tmp1 = vscales(indn+k) * &
5314 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
5315 & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
5318 tmp2 = (two**(-l)) / &
5319 & fact(l+k+1)*fact(2*l+1)/ &
5320 & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)/ &
5324 res1 = res1 + tmp3*src_m(indn+k)
5325 res2 = res2 + tmp3*src_m(indn-k)
5326 dst_m(indj+k) = beta*dst_m(indj+k) + res1
5327 dst_m(indj-k) = beta*dst_m(indj-k) + res2
5340 tmp1 = vscales(indn) * &
5341 & dble(2*j+1) * fact(j+1) * fact(n+1) / &
5342 & src_sk(n+1) / vscales(indj) * src_sk(j+1) * alpha
5345 tmp2 = (two**(-l)) / &
5346 & fact(l+1)*fact(2*l+1)/ &
5347 & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)/ &
5351 res1 = res1 + tmp3*src_m(indn)
5352 dst_m(indj) = beta*dst_m(indj) + res1
5361 tmp1 = vscales(indn+k) * &
5362 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
5363 & src_sk(n+1) / vscales(indj+k) * src_sk(j+1) * alpha
5366 tmp2 = (two**(-l)) / &
5367 & fact(l+k+1)*fact(2*l+1)/ &
5368 & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)/ &
5372 res1 = res1 + tmp3*src_m(indn+k)
5373 res2 = res2 + tmp3*src_m(indn-k)
5374 dst_m(indj+k) = beta*dst_m(indj+k) + res1
5375 dst_m(indj-k) = beta*dst_m(indj-k) + res2
5403 & alpha, src_m, beta, dst_m)
5405 integer,
intent(in) :: p
5406 real(dp),
intent(in) :: z, src_r, dst_r, vscales((p+1)*(p+1)), &
5407 & vcnk((2*p+1)*(p+1)), alpha, src_m((p+1)*(p+1)), beta
5409 real(dp),
intent(inout) :: dst_m((p+1)*(p+1))
5411 real(dp) :: work(2*(p+1))
5414 & alpha, src_m, beta, dst_m, work)
5441 & alpha, src_m, beta, dst_m, work)
5443 integer,
intent(in) :: p
5444 real(dp),
intent(in) :: z, src_r, dst_r, vscales((p+1)*(p+1)), &
5445 & vcnk((2*p+1)*(p+1)), alpha, src_m((p+1)*(p+1)), beta
5447 real(dp),
intent(inout) :: dst_m((p+1)*(p+1))
5449 real(dp),
intent(out),
target :: work(2*(p+1))
5451 real(dp) :: r1, r2, tmp1, tmp2, res1, res2
5452 integer :: j, k, n, indj, indn, indjk1, indjk2
5454 real(dp),
pointer :: pow_r1(:), pow_r2(:)
5456 if (alpha .eq. zero)
then
5457 if (beta .eq. zero)
then
5460 dst_m = beta * dst_m
5467 pow_r1(1:p+1) => work(1:p+1)
5468 pow_r2(1:p+1) => work(p+2:2*(p+1))
5475 pow_r1(j) = pow_r1(j-1) * r1
5476 pow_r2(j) = pow_r2(j-1) * r2
5480 if (beta .eq. zero)
then
5484 tmp1 = alpha * pow_r1(n+1) / vscales(indn)
5491 indjk1 = j*(j+1)/2 + 1
5492 tmp2 = tmp1 * vscales(indj) * pow_r2(j-n+1) * &
5493 & vcnk(indjk1+j-n)**2
5494 res1 = res1 + tmp2*src_m(indj)
5505 indjk1 = (j-k)*(j-k+1)/2 + 1
5506 indjk2 = (j+k)*(j+k+1)/2 + 1
5507 tmp2 = tmp1 * vscales(indj) * pow_r2(j-n+1) * &
5508 & vcnk(indjk1+j-n) * vcnk(indjk2+j-n)
5509 res1 = res1 + tmp2*src_m(indj+k)
5510 res2 = res2 + tmp2*src_m(indj-k)
5512 dst_m(indn+k) = res1
5513 dst_m(indn-k) = res2
5521 tmp1 = alpha * pow_r1(n+1) / vscales(indn)
5528 indjk1 = j*(j+1)/2 + 1
5529 tmp2 = tmp1 * vscales(indj) * pow_r2(j-n+1) * &
5530 & vcnk(indjk1+j-n)**2
5531 res1 = res1 + tmp2*src_m(indj)
5533 dst_m(indn) = beta*dst_m(indn) + res1
5542 indjk1 = (j-k)*(j-k+1)/2 + 1
5543 indjk2 = (j+k)*(j+k+1)/2 + 1
5544 tmp2 = tmp1 * vscales(indj) * pow_r2(j-n+1) * &
5545 & vcnk(indjk1+j-n) * vcnk(indjk2+j-n)
5546 res1 = res1 + tmp2*src_m(indj+k)
5547 res2 = res2 + tmp2*src_m(indj-k)
5549 dst_m(indn+k) = beta*dst_m(indn+k) + res1
5550 dst_m(indn-k) = beta*dst_m(indn-k) + res2
5557 if (beta .eq. zero)
then
5563 do k = indj-j, indj+j
5564 dst_m(k) = tmp1 * src_m(k)
5575 do k = indj-j, indj+j
5576 dst_m(k) = beta*dst_m(k) + tmp1*src_m(k)
5605 integer,
intent(in) :: p
5606 real(dp),
intent(in) :: src_r, dst_r, alpha, src_m((p+1)*(p+1)), beta
5608 real(dp),
intent(inout) :: dst_m((p+1)*(p+1))
5610 real(dp) :: r1, tmp1
5611 integer :: j, k, indj
5613 if (beta .eq. zero)
then
5616 dst_m = beta * dst_m
5622 do k = indj-j, indj+j
5623 dst_m(k) = dst_m(k) + src_m(k)*tmp1
5651 integer,
intent(in) :: p
5652 real(dp),
intent(in) :: src_r, dst_r, alpha, src_m((p+1)*(p+1)), beta
5654 real(dp),
intent(inout) :: dst_m((p+1)*(p+1))
5656 real(dp) :: r1, tmp1
5657 integer :: j, k, indj
5659 if (beta .eq. zero)
then
5662 dst_m = beta * dst_m
5668 do k = indj-j, indj+j
5669 dst_m(k) = dst_m(k) + src_m(k)*tmp1
5701 & src_m, beta, dst_m)
5703 integer,
intent(in) :: p
5704 real(dp),
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
5705 & vcnk((2*p+1)*(p+1)), alpha, src_m((p+1)*(p+1)), beta
5707 real(dp),
intent(inout) :: dst_m((p+1)*(p+1))
5709 real(dp) :: work(6*p*p + 19*p + 8)
5712 & src_m, beta, dst_m, work)
5742 & src_m, beta, dst_m, work)
5744 integer,
intent(in) :: p
5745 real(dp),
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
5746 & vcnk((2*p+1)*(p+1)), alpha, src_m((p+1)*(p+1)), beta
5748 real(dp),
intent(inout) :: dst_m((p+1)*(p+1))
5750 real(dp),
intent(out),
target :: work(6*p*p + 19*p + 8)
5752 real(dp) :: rho, ctheta, stheta, cphi, sphi
5755 real(dp),
pointer :: tmp_m(:), tmp_m2(:), vcos(:), vsin(:)
5757 call carttosph(c, rho, ctheta, stheta, cphi, sphi)
5759 if (stheta .eq. zero)
then
5762 & alpha, src_m, beta, dst_m, work)
5768 tmp_m(1:m) => work(n+1:n+m)
5770 tmp_m2(1:m) => work(n+1:n+m)
5773 vcos => work(n+1:n+m)
5775 vsin => work(n+1:n+m)
5777 call trgev(cphi, sphi, p, vcos, vsin)
5785 & tmp_m2, zero, tmp_m, work)
5818 & src_m, beta, dst_m)
5821 integer,
intent(in) :: p
5822 real(dp),
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
5823 & alpha, src_m((p+1)*(p+1)), beta, kappa
5825 real(dp),
intent(inout) :: dst_m((p+1)*(p+1))
5827 real(dp) :: work(6*p*p + 19*p + 8), src_sk(p+1), dst_sk(p+1), s1, s2, ck(3)
5828 complex(dp) :: work_complex(2*p+1), z1, z2
5833 s1 = sqrt(2/(pi*real(z1)))
5834 s2 = sqrt(2/(pi*real(z2)))
5835 call cbesk(z1, pt5, 1, 1, work_complex(1), nz, ierr)
5836 src_sk(1) = s1 * real(work_complex(1))
5837 call cbesk(z2, pt5, 1, 1, work_complex(1), nz, ierr)
5838 dst_sk(1) = s2 * real(work_complex(1))
5840 call cbesk(z1, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
5841 src_sk(2:p+1) = s1 * real(work_complex(2:p+1))
5842 call cbesk(z2, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
5843 dst_sk(2:p+1) = s2 * real(work_complex(2:p+1))
5848 & src_m, beta, dst_m, work, work_complex)
5877 & src_m, beta, dst_m, work, work_complex)
5879 integer,
intent(in) :: p
5880 real(dp),
intent(in) :: c(3), src_sk(p+1), dst_sk(p+1), vscales((p+1)*(p+1)), &
5881 & alpha, src_m((p+1)*(p+1)), beta
5883 real(dp),
intent(inout) :: dst_m((p+1)*(p+1))
5885 real(dp),
intent(out),
target :: work(6*p*p + 19*p + 8)
5886 complex(dp),
intent(out) :: work_complex(2*p+1)
5888 real(dp) :: rho, ctheta, stheta, cphi, sphi
5891 real(dp),
pointer :: tmp_m(:), tmp_m2(:), vcos(:), vsin(:)
5893 call carttosph(c, rho, ctheta, stheta, cphi, sphi)
5904 tmp_m(1:m) => work(n+1:n+m)
5906 tmp_m2(1:m) => work(n+1:n+m)
5909 vcos => work(n+1:n+m)
5911 vsin => work(n+1:n+m)
5913 call trgev(cphi, sphi, p, vcos, vsin)
5921 & tmp_m2, zero, tmp_m, work, work_complex)
5954 & alpha, src_m, beta, dst_m)
5957 integer,
intent(in) :: p
5958 real(dp),
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
5959 & alpha, src_m((p+1)*(p+1)), beta, kappa
5961 real(dp),
intent(inout) :: dst_m((p+1)*(p+1))
5963 real(dp) :: work(6*p*p + 19*p + 8), src_sk(p+1), dst_sk(p+1), s1, s2, ck(3)
5964 complex(dp) :: work_complex(2*p+1), z1, z2
5969 s1 = sqrt(2/(pi*real(z1)))
5970 s2 = sqrt(2/(pi*real(z2)))
5971 call cbesk(z1, pt5, 1, 1, work_complex(1), nz, ierr)
5972 src_sk(1) = s1 * real(work_complex(1))
5973 call cbesk(z2, pt5, 1, 1, work_complex(1), nz, ierr)
5974 dst_sk(1) = s2 * real(work_complex(1))
5976 call cbesk(z1, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
5977 src_sk(2:p+1) = s1 * real(work_complex(2:p+1))
5978 call cbesk(z2, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
5979 dst_sk(2:p+1) = s2 * real(work_complex(2:p+1))
5984 & alpha, src_m, beta, dst_m, work, work_complex)
6013 & src_m, beta, dst_m, work, work_complex)
6015 integer,
intent(in) :: p
6016 real(dp),
intent(in) :: c(3), src_sk(p+1), dst_sk(p+1), vscales((p+1)*(p+1)), &
6017 & alpha, src_m((p+1)*(p+1)), beta
6019 real(dp),
intent(inout) :: dst_m((p+1)*(p+1))
6021 real(dp),
intent(out),
target :: work(6*p*p + 19*p + 8)
6022 complex(dp),
intent(out) :: work_complex(2*p+1)
6024 real(dp) :: rho, ctheta, stheta, cphi, sphi
6027 real(dp),
pointer :: tmp_m(:), tmp_m2(:), vcos(:), vsin(:)
6029 call carttosph(c, rho, ctheta, stheta, cphi, sphi)
6040 tmp_m(1:m) => work(n+1:n+m)
6042 tmp_m2(1:m) => work(n+1:n+m)
6045 vcos => work(n+1:n+m)
6047 vsin => work(n+1:n+m)
6049 call trgev(cphi, sphi, p, vcos, vsin)
6057 & tmp_m2, zero, tmp_m, work, work_complex)
6089 & src_m, beta, dst_m)
6091 integer,
intent(in) :: p
6092 real(dp),
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
6093 & vcnk((2*p+1)*(p+1)), alpha, src_m((p+1)*(p+1)), beta
6095 real(dp),
intent(inout) :: dst_m((p+1)*(p+1))
6097 real(dp) :: work(6*p*p + 19*p + 8)
6100 & src_m, beta, dst_m, work)
6128 & alpha, src_m, beta, dst_m, work)
6130 integer,
intent(in) :: p
6131 real(dp),
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
6132 & vcnk((2*p+1)*(p+1)), alpha, src_m((p+1)*(p+1)), beta
6134 real(dp),
intent(inout) :: dst_m((p+1)*(p+1))
6136 real(dp),
intent(out),
target :: work(6*p*p + 19*p + 8)
6138 real(dp) :: rho, ctheta, stheta, cphi, sphi
6141 real(dp),
pointer :: tmp_m(:), tmp_m2(:), vcos(:), vsin(:)
6143 call carttosph(c, rho, ctheta, stheta, cphi, sphi)
6145 if (stheta .eq. 0)
then
6148 & vcnk, alpha, src_m, beta, dst_m, work)
6154 tmp_m(1:m) => work(n+1:n+m)
6156 tmp_m2(1:m) => work(n+1:n+m)
6159 vcos => work(n+1:n+m)
6161 vsin => work(n+1:n+m)
6163 call trgev(cphi, sphi, p, vcos, vsin)
6171 & one, tmp_m2, zero, tmp_m, work)
6179subroutine fmm_m2m_bessel_grad(p, sph_sk, vscales, sph_m, sph_m_grad)
6181 integer,
intent(in) :: p
6182 real(dp),
intent(in) :: sph_sk(p+2), vscales((p+2)**2), &
6185 real(dp),
intent(out) :: sph_m_grad((p+2)**2, 3)
6187 real(dp),
dimension(3, 3) :: zx_coord_transform, zy_coord_transform
6188 real(dp) :: tmp((p+1)**2), tmp2((p+2)**2)
6190 zx_coord_transform = zero
6191 zx_coord_transform(3, 2) = one
6192 zx_coord_transform(2, 3) = one
6193 zx_coord_transform(1, 1) = one
6194 zy_coord_transform = zero
6195 zy_coord_transform(1, 2) = one
6196 zy_coord_transform(2, 1) = one
6197 zy_coord_transform(3, 3) = one
6202 & one, tmp, zero, tmp2)
6210 & one, tmp, zero, tmp2)
6216 & one, sph_m, zero, sph_m_grad(:, 3))
6217end subroutine fmm_m2m_bessel_grad
6242 & src_l, beta, dst_l)
6244 integer,
intent(in) :: p
6245 real(dp),
intent(in) :: z, src_r, dst_r, vscales((p+1)*(p+1)), &
6246 & vfact(2*p+1), alpha, src_l((p+1)*(p+1)), beta
6248 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
6250 real(dp) :: work(2*(p+1))
6253 & src_l, beta, dst_l, work)
6280 & src_l, beta, dst_l, work)
6282 integer,
intent(in) :: p
6283 real(dp),
intent(in) :: z, src_r, dst_r, vscales((p+1)*(p+1)), &
6284 & vfact(2*p+1), alpha, src_l((p+1)*(p+1)), beta
6286 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
6288 real(dp),
intent(out),
target :: work(2*(p+1))
6290 real(dp) :: r1, r2, tmp1, tmp2
6291 integer :: j, k, n, indj, indn
6293 real(dp),
pointer :: pow_r1(:), pow_r2(:)
6295 if (beta .eq. zero)
then
6298 dst_l = beta * dst_l
6303 pow_r1(1:p+1) => work(1:p+1)
6304 pow_r2(1:p+1) => work(p+2:2*(p+1))
6311 pow_r1(j) = pow_r1(j-1) * r1
6312 pow_r2(j) = pow_r2(j-1) * r2
6318 tmp1 = alpha * pow_r2(j+1) / vfact(j-k+1) / vfact(j+k+1) * &
6322 tmp2 = tmp1 * pow_r1(n+1) / vscales(indn) * &
6323 & vfact(n-k+1) * vfact(n+k+1) / vfact(n-j+1) / &
6325 if (mod(n+j, 2) .eq. 1)
then
6329 dst_l(indj) = dst_l(indj) + tmp2*src_l(indn)
6331 dst_l(indj+k) = dst_l(indj+k) + tmp2*src_l(indn+k)
6332 dst_l(indj-k) = dst_l(indj-k) + tmp2*src_l(indn-k)
6343 do k = indj-j, indj+j
6344 dst_l(k) = dst_l(k) + src_l(k)*tmp1
6373 & src_l, beta, dst_l)
6375 integer,
intent(in) :: p
6376 real(dp),
intent(in) :: z, src_si(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
6377 & alpha, src_l((p+1)*(p+1)), beta
6379 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
6381 real(dp) :: work(2*p+1)
6382 complex(dp) :: work_complex(2*p+1)
6385 & alpha, src_l, beta, dst_l, work, work_complex)
6412 & alpha, src_l, beta, dst_l, work, work_complex)
6415 integer,
intent(in) :: p
6416 real(dp),
intent(in) :: z, src_si(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
6417 & alpha, src_l((p+1)*(p+1)), beta
6419 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
6421 real(dp),
intent(out),
target :: work(2*p+1)
6422 complex(dp),
intent(out) :: work_complex(2*p+1)
6424 real(dp) :: r1, tmp1, tmp2, tmp3, res1, res2
6425 integer :: j, k, n, indj, indn, l, NZ, ierr
6427 real(dp) :: fact(2*p+1)
6428 complex(dp) :: complex_argument
6430 if (alpha .eq. zero)
then
6431 if (beta .eq. zero)
then
6434 dst_l = beta * dst_l
6441 fact(j+1) = dble(j) * fact(j)
6446 complex_argument = z
6447 call cbesi(complex_argument, pt5, 1, 1, work_complex(1), nz, ierr)
6448 work(1) = sqrt(pi/(2*z)) * real(work_complex(1))
6450 call cbesi(complex_argument, 1.5d0, 1, 2*p, &
6451 & work_complex(2:2*p+1), nz, ierr)
6452 work(2:2*p+1) = sqrt(pi/(2*z)) * real(work_complex(2:2*p+1))
6456 if (beta .eq. zero)
then
6467 tmp1 = vscales(indn) * ((-one)**(j+n)) * &
6468 & dble(2*j+1) * fact(j+1) * fact(n+1) / &
6469 & src_si(n+1) / vscales(indj) * dst_si(j+1) * alpha
6472 tmp2 = (two**(-l)) / &
6473 & fact(l+1)*fact(2*l+1)/ &
6474 & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
6475 & (work(j+n-l+1) / (z**l))
6479 res1 = res1 + tmp3*src_l(indn)
6489 tmp1 = vscales(indn+k) * ((-one)**(j+n)) * &
6490 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
6491 & src_si(n+1) / vscales(indj+k) * dst_si(j+1) * alpha
6494 tmp2 = (two**(-l)) / &
6495 & fact(l+k+1)*fact(2*l+1)/ &
6496 & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
6497 & (work(j+n-l+1) / (z**l))
6501 res1 = res1 + tmp3*src_l(indn+k)
6502 res2 = res2 + tmp3*src_l(indn-k)
6504 dst_l(indj+k) = res1
6505 dst_l(indj-k) = res2
6519 tmp1 = vscales(indn) * ((-one)**(j+n)) * &
6520 & dble(2*j+1) * fact(j+1) * fact(n+1) / &
6521 & src_si(n+1) / vscales(indj) * dst_si(j+1) * alpha
6524 tmp2 = (two**(-l)) / &
6525 & fact(l+1)*fact(2*l+1)/ &
6526 & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
6527 & (work(j+n-l+1) / (z**l))
6531 res1 = res1 + tmp3*src_l(indn)
6533 dst_l(indj) = beta*dst_l(indj) + res1
6541 tmp1 = vscales(indn+k) * ((-one)**(j+n)) * &
6542 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
6543 & src_si(n+1) / vscales(indj+k) * dst_si(j+1) * alpha
6546 tmp2 = (two**(-l)) / &
6547 & fact(l+k+1)*fact(2*l+1)/ &
6548 & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
6549 & (work(j+n-l+1) / (z**l))
6553 res1 = res1 + tmp3*src_l(indn+k)
6554 res2 = res2 + tmp3*src_l(indn-k)
6556 dst_l(indj+k) = beta*dst_l(indj+k) + res1
6557 dst_l(indj-k) = beta*dst_l(indj-k) + res2
6565 if (beta .eq. zero)
then
6570 do k = indj-j, indj+j
6571 dst_l(k) = tmp1 * src_l(k)
6581 do k = indj-j, indj+j
6582 dst_l(k) = beta*dst_l(k) + tmp1*src_l(k)
6612 & src_l, beta, dst_l)
6614 integer,
intent(in) :: p
6615 real(dp),
intent(in) :: z, src_si(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
6616 & alpha, src_l((p+1)*(p+1)), beta
6618 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
6620 real(dp) :: work(2*p+1)
6621 complex(dp) :: work_complex(2*p+1)
6624 & alpha, src_l, beta, dst_l, work, work_complex)
6651 & alpha, src_l, beta, dst_l, work, work_complex)
6654 integer,
intent(in) :: p
6655 real(dp),
intent(in) :: z, src_si(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
6656 & alpha, src_l((p+1)*(p+1)), beta
6658 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
6660 real(dp),
intent(out),
target :: work(2*p+1)
6661 complex(dp),
intent(out) :: work_complex(2*p+1)
6663 real(dp) :: r1, tmp1, tmp2, tmp3, res1, res2
6664 integer :: j, k, n, indj, indn, l, NZ, ierr
6666 real(dp) :: fact(2*p+1)
6667 complex(dp) :: complex_argument
6669 if (alpha .eq. zero)
then
6670 if (beta .eq. zero)
then
6673 dst_l = beta * dst_l
6680 fact(j+1) = dble(j) * fact(j)
6685 complex_argument = z
6686 call cbesi(complex_argument, pt5, 1, 1, work_complex(1), nz, ierr)
6687 work(1) = sqrt(pi/(2*z)) * real(work_complex(1))
6689 call cbesi(complex_argument, 1.5d0, 1, 2*p, &
6690 & work_complex(2:2*p+1), nz, ierr)
6691 work(2:2*p+1) = sqrt(pi/(2*z)) * real(work_complex(2:2*p+1))
6695 if (beta .eq. zero)
then
6706 tmp1 = vscales(indn) * ((-one)**(j+n)) * &
6707 & dble(2*j+1) * fact(j+1) * fact(n+1) / &
6708 & src_si(n+1) / vscales(indj) * dst_si(j+1) * alpha
6711 tmp2 = (two**(-l)) / &
6712 & fact(l+1)*fact(2*l+1)/ &
6713 & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
6714 & (work(j+n-l+1) / (z**l))
6718 res1 = res1 + tmp3*src_l(indj)
6728 tmp1 = vscales(indn+k) * ((-one)**(j+n)) * &
6729 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
6730 & src_si(n+1) / vscales(indj+k) * dst_si(j+1) * alpha
6733 tmp2 = (two**(-l)) / &
6734 & fact(l+k+1)*fact(2*l+1)/ &
6735 & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
6736 & (work(j+n-l+1) / (z**l))
6740 res1 = res1 + tmp3*src_l(indj+k)
6741 res2 = res2 + tmp3*src_l(indj-k)
6743 dst_l(indn+k) = res1
6744 dst_l(indn-k) = res2
6758 tmp1 = vscales(indn) * ((-one)**(j+n)) * &
6759 & dble(2*j+1) * fact(j+1) * fact(n+1) / &
6761 & src_si(n+1) / vscales(indj) * dst_si(j+1) * alpha
6764 tmp2 = (two**(-l)) / &
6765 & fact(l+1)*fact(2*l+1)/ &
6766 & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
6767 & (work(j+n-l+1) / (z**l))
6771 res1 = res1 + tmp3*src_l(indj)
6773 dst_l(indn) = beta*dst_l(indn) + res1
6781 tmp1 = vscales(indn+k) * ((-one)**(j+n)) * &
6782 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
6783 & src_si(n+1) / vscales(indj+k) * dst_si(j+1) * alpha
6786 tmp2 = (two**(-l)) / &
6787 & fact(l+k+1)*fact(2*l+1)/ &
6788 & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
6789 & (work(j+n-l+1) / (z**l))
6793 res1 = res1 + tmp3*src_l(indj+k)
6794 res2 = res2 + tmp3*src_l(indj-k)
6796 dst_l(indn+k) = beta*dst_l(indn+k) + res1
6797 dst_l(indn-k) = beta*dst_l(indn-k) + res2
6805 if (beta .eq. zero)
then
6810 do k = indj-j, indj+j
6811 dst_l(k) = tmp1 * src_l(k)
6821 do k = indj-j, indj+j
6822 dst_l(k) = beta*dst_l(k) + tmp1*src_l(k)
6854 & alpha, src_l, beta, dst_l)
6857 integer,
intent(in) :: p
6858 real(dp),
intent(in) :: src_si(p+2), vscales((p+2)*(p+2)), &
6859 & alpha, src_l((p+1)*(p+1)), beta
6861 real(dp),
intent(inout) :: dst_l((p+2)*(p+2))
6864 & -alpha, src_l, beta, dst_l)
6890 & alpha, src_l, beta, dst_l)
6892 integer,
intent(in) :: p
6893 real(dp),
intent(in) :: z, src_r, dst_r, vscales((p+1)*(p+1)), &
6894 & vfact(2*p+1), alpha, src_l((p+1)*(p+1)), beta
6896 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
6898 real(dp) :: work(2*(p+1))
6901 & alpha, src_l, beta, dst_l, work)
6928 & alpha, src_l, beta, dst_l, work)
6930 integer,
intent(in) :: p
6931 real(dp),
intent(in) :: z, src_r, dst_r, vscales((p+1)*(p+1)), &
6932 & vfact(2*p+1), alpha, src_l((p+1)*(p+1)), beta
6934 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
6936 real(dp),
intent(out),
target :: work(2*(p+1))
6938 real(dp) :: r1, r2, tmp1, tmp2
6939 integer :: j, k, n, indj, indn
6941 real(dp),
pointer :: pow_r1(:), pow_r2(:)
6943 if (beta .eq. zero)
then
6946 dst_l = beta * dst_l
6951 pow_r1(1:p+1) => work(1:p+1)
6952 pow_r2(1:p+1) => work(p+2:2*(p+1))
6959 pow_r1(j) = pow_r1(j-1) * r1
6960 pow_r2(j) = pow_r2(j-1) * r2
6965 tmp1 = alpha * pow_r2(j+1) / vfact(j-k+1) / vfact(j+k+1) * &
6969 tmp2 = tmp1 * pow_r1(n+1) / vscales(indn) * &
6970 & vfact(n-k+1) * vfact(n+k+1) / vfact(n-j+1) / &
6972 if (mod(n+j, 2) .eq. 1)
then
6976 dst_l(indn) = dst_l(indn) + tmp2*src_l(indj)
6978 dst_l(indn+k) = dst_l(indn+k) + tmp2*src_l(indj+k)
6979 dst_l(indn-k) = dst_l(indn-k) + tmp2*src_l(indj-k)
6990 do k = indj-j, indj+j
6991 dst_l(k) = dst_l(k) + src_l(k)*tmp1
7019 integer,
intent(in) :: p
7020 real(dp),
intent(in) :: src_r, dst_r, alpha, src_l((p+1)*(p+1)), beta
7022 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
7024 real(dp) :: r1, tmp1
7025 integer :: j, k, indj
7027 if (beta .eq. zero)
then
7030 dst_l = beta * dst_l
7036 do k = indj-j, indj+j
7037 dst_l(k) = dst_l(k) + src_l(k)*tmp1
7064 integer,
intent(in) :: p
7065 real(dp),
intent(in) :: src_r, dst_r, alpha, src_l((p+1)*(p+1)), beta
7067 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
7069 real(dp) :: r1, tmp1
7070 integer :: j, k, indj
7072 if (beta .eq. zero)
then
7075 dst_l = beta * dst_l
7081 do k = indj-j, indj+j
7082 dst_l(k) = dst_l(k) + src_l(k)*tmp1
7114 & src_l, beta, dst_l)
7116 integer,
intent(in) :: p
7117 real(dp),
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
7118 & vfact(2*p+1), alpha, src_l((p+1)*(p+1)), beta
7120 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
7122 real(dp) :: work(6*p*p+19*p+8)
7125 & src_l, beta, dst_l, work)
7155 & src_l, beta, dst_l, work)
7157 integer,
intent(in) :: p
7158 real(dp),
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
7159 & vfact(2*p+1), alpha, src_l((p+1)*(p+1)), beta
7161 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
7163 real(dp),
intent(out),
target :: work(6*p*p + 19*p + 8)
7165 real(dp) :: rho, ctheta, stheta, cphi, sphi
7168 real(dp),
pointer :: tmp_l(:), tmp_l2(:), vcos(:), vsin(:)
7170 call carttosph(c, rho, ctheta, stheta, cphi, sphi)
7172 if (stheta .eq. zero)
then
7175 & alpha, src_l, beta, dst_l, work)
7181 tmp_l(1:m) => work(n+1:n+m)
7183 tmp_l2(1:m) => work(n+1:n+m)
7186 vcos => work(n+1:n+m)
7188 vsin => work(n+1:n+m)
7190 call trgev(cphi, sphi, p, vcos, vsin)
7198 & tmp_l2, zero, tmp_l, work)
7231 & src_l, beta, dst_l)
7234 integer,
intent(in) :: p
7235 real(dp),
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
7236 & alpha, src_l((p+1)*(p+1)), beta, kappa
7238 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
7240 real(dp) :: work(6*p*p + 19*p + 8), src_si(p+1), dst_si(p+1), s1, s2
7241 complex(dp) :: work_complex(2*p+1), z1, z2
7246 s1 = sqrt(pi/(2*real(z1)))
7247 s2 = sqrt(pi/(2*real(z2)))
7248 call cbesi(z1, pt5, 1, 1, work_complex(1), nz, ierr)
7249 src_si(1) = s1 * real(work_complex(1))
7250 call cbesi(z2, pt5, 1, 1, work_complex(1), nz, ierr)
7251 dst_si(1) = s2 * real(work_complex(1))
7253 call cbesi(z1, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
7254 src_si(2:p+1) = s1 * real(work_complex(2:p+1))
7255 call cbesi(z2, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
7256 dst_si(2:p+1) = s2 * real(work_complex(2:p+1))
7260 & src_l, beta, dst_l, work, work_complex)
7289 & src_l, beta, dst_l, work, work_complex)
7291 integer,
intent(in) :: p
7292 real(dp),
intent(in) :: c(3), src_si(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
7293 & alpha, src_l((p+1)*(p+1)), beta
7295 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
7297 real(dp),
intent(out),
target :: work(6*p*p + 19*p + 8)
7298 complex(dp),
intent(out) :: work_complex(2*p+1)
7300 real(dp) :: rho, ctheta, stheta, cphi, sphi
7303 real(dp),
pointer :: tmp_l(:), tmp_l2(:), vcos(:), vsin(:)
7305 call carttosph(c, rho, ctheta, stheta, cphi, sphi)
7315 tmp_l(1:m) => work(n+1:n+m)
7317 tmp_l2(1:m) => work(n+1:n+m)
7320 vcos => work(n+1:n+m)
7322 vsin => work(n+1:n+m)
7324 call trgev(cphi, sphi, p, vcos, vsin)
7332 & one, tmp_l2, zero, tmp_l, work, work_complex)
7365 & src_l, beta, dst_l)
7368 integer,
intent(in) :: p
7369 real(dp),
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
7370 & alpha, src_l((p+1)*(p+1)), beta, kappa
7372 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
7374 real(dp) :: work(6*p*p + 19*p + 8), src_si(p+1), dst_si(p+1), s1, s2
7375 complex(dp) :: work_complex(2*p+1), z1, z2
7380 s1 = sqrt(pi/(2*real(z1)))
7381 s2 = sqrt(pi/(2*real(z2)))
7382 call cbesi(z1, pt5, 1, 1, work_complex(1), nz, ierr)
7383 src_si(1) = s1 * real(work_complex(1))
7384 call cbesi(z2, pt5, 1, 1, work_complex(1), nz, ierr)
7385 dst_si(1) = s2 * real(work_complex(1))
7387 call cbesi(z1, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
7388 src_si(2:p+1) = s1 * real(work_complex(2:p+1))
7389 call cbesi(z2, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
7390 dst_si(2:p+1) = s2 * real(work_complex(2:p+1))
7394 & src_l, beta, dst_l, work, work_complex)
7423 & alpha, src_l, beta, dst_l, work, work_complex)
7425 integer,
intent(in) :: p
7426 real(dp),
intent(in) :: c(3), src_si(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
7427 & alpha, src_l((p+1)*(p+1)), beta
7429 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
7431 real(dp),
intent(out),
target :: work(6*p*p + 19*p + 8)
7432 complex(dp),
intent(out) :: work_complex(2*p+1)
7434 real(dp) :: rho, ctheta, stheta, cphi, sphi
7437 real(dp),
pointer :: tmp_l(:), tmp_l2(:), vcos(:), vsin(:)
7439 call carttosph(c, rho, ctheta, stheta, cphi, sphi)
7450 tmp_l(1:m) => work(n+1:n+m)
7452 tmp_l2(1:m) => work(n+1:n+m)
7455 vcos => work(n+1:n+m)
7457 vsin => work(n+1:n+m)
7459 call trgev(cphi, sphi, p, vcos, vsin)
7467 & one, tmp_l2, zero, tmp_l, work, work_complex)
7475subroutine fmm_l2l_bessel_grad(p, sph_si, vscales, sph_l, sph_l_grad)
7477 integer,
intent(in) :: p
7478 real(dp),
intent(in) :: sph_si(p+2), vscales((p+2)**2), sph_l((p+1)**2)
7480 real(dp),
intent(out) :: sph_l_grad((p+2)**2, 3)
7482 real(dp),
dimension(3, 3) :: zx_coord_transform, zy_coord_transform
7483 real(dp) :: tmp((p+1)**2), tmp2((p+2)**2)
7485 zx_coord_transform = zero
7486 zx_coord_transform(3, 2) = one
7487 zx_coord_transform(2, 3) = one
7488 zx_coord_transform(1, 1) = one
7489 zy_coord_transform = zero
7490 zy_coord_transform(1, 2) = one
7491 zy_coord_transform(2, 1) = one
7492 zy_coord_transform(3, 3) = one
7497 & one, tmp, zero, tmp2)
7505 & one, tmp, zero, tmp2)
7511 & one, sph_l, zero, sph_l_grad(:, 3))
7512end subroutine fmm_l2l_bessel_grad
7540 & alpha, src_l, beta, dst_l)
7542 integer,
intent(in) :: p
7543 real(dp),
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
7544 & vfact(2*p+1), alpha, src_l((p+1)*(p+1)), beta
7546 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
7548 real(dp) :: work(6*p*p + 19*p + 8)
7551 & src_l, beta, dst_l, work)
7581 & alpha, src_l, beta, dst_l, work)
7583 integer,
intent(in) :: p
7584 real(dp),
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
7585 & vfact(2*p+1), alpha, src_l((p+1)*(p+1)), beta
7587 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
7589 real(dp),
intent(out),
target :: work(6*p*p + 19*p + 8)
7591 real(dp) :: rho, ctheta, stheta, cphi, sphi
7594 real(dp),
pointer :: tmp_l(:), tmp_l2(:), vcos(:), vsin(:)
7596 call carttosph(c, rho, ctheta, stheta, cphi, sphi)
7598 if (stheta .eq. zero)
then
7601 & vfact, alpha, src_l, beta, dst_l, work)
7607 tmp_l(1:m) => work(n+1:n+m)
7609 tmp_l2(1:m) => work(n+1:n+m)
7612 vcos => work(n+1:n+m)
7614 vsin => work(n+1:n+m)
7616 call trgev(cphi, sphi, p, vcos, vsin)
7624 & one, tmp_l2, zero, tmp_l, work)
7656 & m2l_ztranslate_coef, alpha, src_m, beta, dst_l)
7658 integer,
intent(in) :: pm, pl
7659 real(dp),
intent(in) :: z, src_r, dst_r, vscales((pm+pl+1)*(pm+pl+1)), &
7660 & m2l_ztranslate_coef(pm+1, pl+1, pl+1), alpha, src_m((pm+1)*(pm+1)), &
7663 real(dp),
intent(inout) :: dst_l((pl+1)*(pl+1))
7665 real(dp) :: work((pm+2)*(pm+1))
7668 & m2l_ztranslate_coef, alpha, src_m, beta, dst_l, work)
7696 & m2l_ztranslate_coef, alpha, src_m, beta, dst_l, work)
7698 integer,
intent(in) :: pm, pl
7699 real(dp),
intent(in) :: z, src_r, dst_r, vscales((pm+pl+1)*(pm+pl+1)), &
7700 & m2l_ztranslate_coef(pm+1, pl+1, pl+1), alpha, src_m((pm+1)*(pm+1)), &
7703 real(dp),
intent(inout) :: dst_l((pl+1)*(pl+1))
7705 real(dp),
intent(out),
target :: work((pm+2)*(pm+1))
7707 real(dp) :: tmp1, r1, r2, res1, res2, pow_r2
7708 integer :: j, k, n, indj, indk1, indk2
7710 real(dp),
pointer :: src_m2(:), pow_r1(:)
7712 if (alpha .eq. zero)
then
7713 if (beta .eq. zero)
then
7716 dst_l = beta * dst_l
7723 if (z .eq. zero)
then
7728 src_m2(1:n) => work(1:n)
7729 pow_r1(1:pm+1) => work(n+1:n+pm+1)
7738 pow_r1(j+1) = tmp1 / vscales(indj)
7746 src_m2(j+1) = pow_r1(j+1) * src_m(indj)
7756 src_m2(indk1+j-k) = pow_r1(j+1) * src_m(indj+k)
7757 src_m2(indk2+j-k) = pow_r1(j+1) * src_m(indj-k)
7763 if (beta .eq. zero)
then
7768 tmp1 = alpha * vscales(indj) * pow_r2
7769 pow_r2 = pow_r2 * r2
7773 res1 = res1 + m2l_ztranslate_coef(n+1, 1, j+1)*src_m2(n+1)
7775 dst_l(indj) = tmp1 * res1
7780 indk1 = pm + 2 + (2*pm-k+2)*(k-1)
7781 indk2 = indk1 + pm - k + 1
7787 & m2l_ztranslate_coef(n-k+1, k+1, j-k+1)* &
7790 & m2l_ztranslate_coef(n-k+1, k+1, j-k+1)* &
7793 dst_l(indj+k) = tmp1 * res1
7794 dst_l(indj-k) = tmp1 * res2
7802 tmp1 = alpha * vscales(indj) * pow_r2
7803 pow_r2 = pow_r2 * r2
7806 res1 = res1 + m2l_ztranslate_coef(n+1, 1, j+1)*src_m2(n+1)
7808 dst_l(indj) = beta*dst_l(indj) + tmp1*res1
7812 indk1 = pm + 2 + (2*pm-k+2)*(k-1)
7813 indk2 = indk1 + pm - k + 1
7818 & m2l_ztranslate_coef(n-k+1, k+1, j-k+1)* &
7821 & m2l_ztranslate_coef(n-k+1, k+1, j-k+1)* &
7824 dst_l(indj+k) = beta*dst_l(indj+k) + tmp1*res1
7825 dst_l(indj-k) = beta*dst_l(indj-k) + tmp1*res2
7855 & m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m)
7857 integer,
intent(in) :: pl, pm
7858 real(dp),
intent(in) :: z, src_r, dst_r, vscales((pm+pl+1)*(pm+pl+1)), &
7859 & m2l_ztranslate_adj_coef(pl+1, pl+1, pm+1), alpha, &
7860 & src_l((pl+1)*(pl+1)), beta
7862 real(dp),
intent(inout) :: dst_m((pm+1)*(pm+1))
7864 real(dp) :: work((pl+2)*(pl+1))
7867 & m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m, work)
7894 & m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m, work)
7896 integer,
intent(in) :: pl, pm
7897 real(dp),
intent(in) :: z, src_r, dst_r, vscales((pm+pl+1)*(pm+pl+1)), &
7898 & m2l_ztranslate_adj_coef(pl+1, pl+1, pm+1), alpha, &
7899 & src_l((pl+1)*(pl+1)), beta
7901 real(dp),
intent(inout) :: dst_m((pm+1)*(pm+1))
7903 real(dp),
intent(out),
target :: work((pl+2)*(pl+1))
7905 real(dp) :: tmp1, r1, r2, pow_r1, res1, res2
7906 integer :: j, k, n, indj, indn, indk1, indk2
7908 real(dp),
pointer :: pow_r2(:), src_l2(:)
7910 if (alpha .eq. zero)
then
7911 if (beta .eq. zero)
then
7914 dst_m = beta * dst_m
7921 if (z .eq. zero)
then
7926 src_l2(1:n) => work(1:n)
7927 pow_r2(1:pl+1) => work(n+1:n+pl+1)
7936 pow_r2(j+1) = tmp1 * vscales(indj)
7944 src_l2(j+1) = pow_r2(j+1) * src_l(indj)
7953 src_l2(indk1+j-k) = pow_r2(j+1) * src_l(indj+k)
7954 src_l2(indk2+j-k) = pow_r2(j+1) * src_l(indj-k)
7960 if (beta .eq. zero)
then
7964 tmp1 = alpha * pow_r1 / vscales(indn)
7965 pow_r1 = pow_r1 * r1
7969 res1 = res1 + m2l_ztranslate_adj_coef(j+1, 1, n+1)*src_l2(j+1)
7971 dst_m(indn) = tmp1 * res1
7975 indk1 = pl + 2 + (2*pl-k+2)*(k-1)
7976 indk2 = indk1 + pl - k + 1
7982 & m2l_ztranslate_adj_coef(j-k+1, k+1, n-k+1)* &
7985 & m2l_ztranslate_adj_coef(j-k+1, k+1, n-k+1)* &
7988 dst_m(indn+k) = tmp1 * res1
7989 dst_m(indn-k) = tmp1 * res2
7997 tmp1 = alpha * pow_r1 / vscales(indn)
7998 pow_r1 = pow_r1 * r1
8002 res1 = res1 + m2l_ztranslate_adj_coef(j+1, 1, n+1)*src_l2(j+1)
8004 dst_m(indn) = beta*dst_m(indn) + tmp1*res1
8008 indk1 = pl + 2 + (2*pl-k+2)*(k-1)
8009 indk2 = indk1 + pl - k + 1
8015 & m2l_ztranslate_adj_coef(j-k+1, k+1, n-k+1)* &
8018 & m2l_ztranslate_adj_coef(j-k+1, k+1, n-k+1)* &
8021 dst_m(indn+k) = beta*dst_m(indn+k) + tmp1*res1
8022 dst_m(indn-k) = beta*dst_m(indn-k) + tmp1*res2
8055 & m2l_ztranslate_coef, alpha, src_m, beta, dst_l)
8057 integer,
intent(in) :: pm, pl
8058 real(dp),
intent(in) :: c(3), src_r, dst_r, vscales((pm+pl+1)**2), &
8059 & m2l_ztranslate_coef(pm+1, pl+1, pl+1), alpha, src_m((pm+1)*(pm+1)), &
8062 real(dp),
intent(inout) :: dst_l((pl+1)*(pl+1))
8064 real(dp) :: work(6*max(pm, pl)**2 + 19*max(pm, pl) + 8)
8067 & m2l_ztranslate_coef, alpha, src_m, beta, dst_l, work)
8099 & m2l_ztranslate_coef, alpha, src_m, beta, dst_l, work)
8101 integer,
intent(in) :: pm, pl
8102 real(dp),
intent(in) :: c(3), src_r, dst_r, vscales((pm+pl+1)**2), &
8103 & m2l_ztranslate_coef(pm+1, pl+1, pl+1), alpha, src_m((pm+1)*(pm+1)), &
8106 real(dp),
intent(inout) :: dst_l((pl+1)*(pl+1))
8108 real(dp),
intent(out),
target :: &
8109 & work(6*max(pm, pl)**2 + 19*max(pm, pl) + 8)
8111 real(dp) :: rho, ctheta, stheta, cphi, sphi
8114 real(dp),
pointer :: tmp_ml(:), tmp_ml2(:), vcos(:), vsin(:)
8116 call carttosph(c, rho, ctheta, stheta, cphi, sphi)
8118 if (stheta .eq. zero)
then
8121 & m2l_ztranslate_coef, alpha, src_m, beta, dst_l, work)
8128 tmp_ml(1:m) => work(n+1:n+m)
8130 tmp_ml2(1:m) => work(n+1:n+m)
8133 vcos => work(n+1:n+m)
8135 vsin => work(n+1:n+m)
8137 call trgev(cphi, sphi, p, vcos, vsin)
8145 & m2l_ztranslate_coef, one, tmp_ml2, zero, tmp_ml, work)
8178 & src_m, beta, dst_l)
8181 integer,
intent(in) :: p
8182 real(dp),
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
8183 & alpha, src_m((p+1)*(p+1)), beta, kappa
8185 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
8187 real(dp) :: work(6*p*p + 19*p + 8), src_sk(p+1), dst_si(p+1), s1, s2
8188 complex(dp) :: work_complex(2*p+1), z1, z2
8193 s1 = sqrt(2/(pi*real(z1)))
8194 s2 = sqrt(pi/(2*real(z2)))
8195 call cbesk(z1, pt5, 1, 1, work_complex(1), nz, ierr)
8196 src_sk(1) = s1 * real(work_complex(1))
8197 call cbesi(z2, pt5, 1, 1, work_complex(1), nz, ierr)
8198 dst_si(1) = s2 * real(work_complex(1))
8200 call cbesk(z1, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
8201 src_sk(2:p+1) = s1 * real(work_complex(2:p+1))
8202 call cbesi(z2, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
8203 dst_si(2:p+1) = s2 * real(work_complex(2:p+1))
8207 & src_m, beta, dst_l, work, work_complex)
8236 & src_m, beta, dst_l, work, work_complex)
8238 integer,
intent(in) :: p
8239 real(dp),
intent(in) :: c(3), src_sk(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
8240 & alpha, src_m((p+1)*(p+1)), beta
8242 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
8244 real(dp),
intent(out),
target :: work(6*p*p + 19*p + 8)
8245 complex(dp),
intent(out) :: work_complex(2*p+1)
8247 real(dp) :: rho, ctheta, stheta, cphi, sphi
8250 real(dp),
pointer :: tmp_m(:), tmp_m2(:), vcos(:), vsin(:)
8252 call carttosph(c, rho, ctheta, stheta, cphi, sphi)
8263 tmp_m(1:m) => work(n+1:n+m)
8265 tmp_m2(1:m) => work(n+1:n+m)
8268 vcos => work(n+1:n+m)
8270 vsin => work(n+1:n+m)
8272 call trgev(cphi, sphi, p, vcos, vsin)
8280 & tmp_m2, zero, tmp_m, work, work_complex)
8313 & src_m, beta, dst_l)
8316 integer,
intent(in) :: p
8317 real(dp),
intent(in) :: c(3), src_r, dst_r, vscales((p+1)*(p+1)), &
8318 & alpha, src_m((p+1)*(p+1)), beta, kappa
8320 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
8322 real(dp) :: work(6*p*p + 19*p + 8), src_sk(p+1), dst_si(p+1), s1, s2, ck(3)
8323 complex(dp) :: work_complex(2*p+1), z1, z2
8328 s1 = sqrt(2/(pi*real(z1)))
8329 s2 = sqrt(pi/(2*real(z2)))
8330 call cbesk(z1, pt5, 1, 1, work_complex(1), nz, ierr)
8331 src_sk(1) = s1 * real(work_complex(1))
8332 call cbesi(z2, pt5, 1, 1, work_complex(1), nz, ierr)
8333 dst_si(1) = s2 * real(work_complex(1))
8335 call cbesk(z1, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
8336 src_sk(2:p+1) = s1 * real(work_complex(2:p+1))
8337 call cbesi(z2, 1.5d0, 1, p, work_complex(2:p+1), nz, ierr)
8338 dst_si(2:p+1) = s2 * real(work_complex(2:p+1))
8343 & alpha, src_m, beta, dst_l, work, work_complex)
8372 & src_m, beta, dst_l, work, work_complex)
8374 integer,
intent(in) :: p
8375 real(dp),
intent(in) :: c(3), src_sk(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
8376 & alpha, src_m((p+1)*(p+1)), beta
8378 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
8380 real(dp),
intent(out),
target :: work(6*p*p + 19*p + 8)
8381 complex(dp),
intent(out) :: work_complex(2*p+1)
8383 real(dp) :: rho, ctheta, stheta, cphi, sphi
8386 real(dp),
pointer :: tmp_m(:), tmp_m2(:), vcos(:), vsin(:)
8388 call carttosph(c, rho, ctheta, stheta, cphi, sphi)
8399 tmp_m(1:m) => work(n+1:n+m)
8401 tmp_m2(1:m) => work(n+1:n+m)
8404 vcos => work(n+1:n+m)
8406 vsin => work(n+1:n+m)
8408 call trgev(cphi, sphi, p, vcos, vsin)
8416 & tmp_m2, zero, tmp_m, work, work_complex)
8446 & src_m, beta, dst_l)
8448 integer,
intent(in) :: p
8449 real(dp),
intent(in) :: z, src_sk(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
8450 & alpha, src_m((p+1)*(p+1)), beta
8452 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
8454 real(dp) :: work(2*p+1)
8455 complex(dp) :: work_complex(2*p+1)
8458 & alpha, src_m, beta, dst_l, work, work_complex)
8485 & alpha, src_m, beta, dst_l, work, work_complex)
8488 integer,
intent(in) :: p
8489 real(dp),
intent(in) :: z, src_sk(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
8490 & alpha, src_m((p+1)*(p+1)), beta
8492 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
8494 real(dp),
intent(out),
target :: work(2*p+1)
8495 complex(dp),
intent(out) :: work_complex(2*p+1)
8497 real(dp) :: r1, tmp1, tmp2, tmp3, res1, res2
8498 integer :: j, k, n, indj, indn, l, NZ, ierr
8500 real(dp) :: fact(2*p+1)
8501 complex(dp) :: complex_argument
8503 if (alpha .eq. zero)
then
8504 if (beta .eq. zero)
then
8507 dst_l = beta * dst_l
8514 fact(j+1) = dble(j) * fact(j)
8519 complex_argument = z
8520 call cbesk(complex_argument, pt5, 1, 1, work_complex(1), nz, ierr)
8521 work(1) = sqrt(2/(pi*z)) * real(work_complex(1))
8523 call cbesk(complex_argument, 1.5d0, 1, 2*p, &
8524 & work_complex(2:2*p+1), nz, ierr)
8525 work(2:2*p+1) = sqrt(2/(pi*z)) * real(work_complex(2:2*p+1))
8529 if (beta .eq. zero)
then
8540 tmp1 = vscales(indn) * ((-one)**(n)) * &
8541 & dble(2*j+1) * fact(j+1) * fact(n+1) / &
8542 & src_sk(n+1) / vscales(indj) * dst_si(j+1)
8545 tmp2 = ((-two)**(-l)) / &
8546 & fact(l+1)*fact(2*l+1)/ &
8547 & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
8548 & (work(j+n-l+1) / (z**l))
8552 res1 = res1 + tmp3*src_m(indn)
8562 tmp1 = vscales(indn+k) * ((-one)**(n)) * &
8563 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
8564 & src_sk(n+1) / vscales(indj+k) * dst_si(j+1)
8567 tmp2 = ((-two)**(-l)) / &
8568 & fact(l+k+1)*fact(2*l+1)/ &
8569 & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
8570 & (work(j+n-l+1) / (z**l))
8574 res1 = res1 + tmp3*src_m(indn+k)
8575 res2 = res2 + tmp3*src_m(indn-k)
8577 dst_l(indj+k) = res1
8578 dst_l(indj-k) = res2
8592 tmp1 = vscales(indn) * ((-one)**(n)) * &
8593 & dble(2*j+1) * fact(j+1) * fact(n+1) / &
8594 & src_sk(n+1) / vscales(indj) * dst_si(j+1)
8597 tmp2 = ((-two)**(-l)) / &
8598 & fact(l+1)*fact(2*l+1)/ &
8599 & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
8600 & (work(j+n-l+1) / (z**l))
8604 res1 = res1 + tmp3*src_m(indn)
8606 dst_l(indj) = beta*dst_l(indj) + res1
8614 tmp1 = vscales(indn+k) * ((-one)**(n)) * &
8615 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) / &
8616 & src_sk(n+1) / vscales(indj+k) * dst_si(j+1)
8619 tmp2 = ((-two)**(-l)) / &
8620 & fact(l+k+1)*fact(2*l+1)/ &
8621 & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
8622 & (work(j+n-l+1) / (z**l))
8626 res1 = res1 + tmp3*src_m(indn+k)
8627 res2 = res2 + tmp3*src_m(indn-k)
8629 dst_l(indj+k) = beta*dst_l(indj+k) + res1
8630 dst_l(indj-k) = beta*dst_l(indj-k) + res2
8638 if (beta .eq. zero)
then
8643 do k = indj-j, indj+j
8644 dst_l(k) = tmp1 * src_m(k)
8654 do k = indj-j, indj+j
8655 dst_l(k) = beta*dst_l(k) + tmp1*src_m(k)
8685 & src_m, beta, dst_l)
8687 integer,
intent(in) :: p
8688 real(dp),
intent(in) :: z, src_sk(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
8689 & alpha, src_m((p+1)*(p+1)), beta
8691 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
8693 real(dp) :: work(2*p+1)
8694 complex(dp) :: work_complex(2*p+1)
8697 & alpha, src_m, beta, dst_l, work, work_complex)
8724 & alpha, src_m, beta, dst_l, work, work_complex)
8727 integer,
intent(in) :: p
8728 real(dp),
intent(in) :: z, src_sk(p+1), dst_si(p+1), vscales((p+1)*(p+1)), &
8729 & alpha, src_m((p+1)*(p+1)), beta
8731 real(dp),
intent(inout) :: dst_l((p+1)*(p+1))
8733 real(dp),
intent(out),
target :: work(2*p+1)
8734 complex(dp),
intent(out) :: work_complex(2*p+1)
8736 real(dp) :: r1, tmp1, tmp2, tmp3, res1, res2
8737 integer :: j, k, n, indj, indn, l, NZ, ierr
8739 real(dp) :: fact(2*p+1)
8740 complex(dp) :: complex_argument
8742 if (alpha .eq. zero)
then
8743 if (beta .eq. zero)
then
8746 dst_l = beta * dst_l
8753 fact(j+1) = dble(j) * fact(j)
8758 complex_argument = z
8759 call cbesk(complex_argument, pt5, 1, 1, work_complex(1), nz, ierr)
8760 work(1) = sqrt(2/(pi*z)) * real(work_complex(1))
8762 call cbesk(complex_argument, 1.5d0, 1, 2*p, &
8763 & work_complex(2:2*p+1), nz, ierr)
8764 work(2:2*p+1) = sqrt(2/(pi*z)) * real(work_complex(2:2*p+1))
8768 if (beta .eq. zero)
then
8779 tmp1 = vscales(indn) * ((-one)**(n)) * &
8780 & dble(2*j+1) * fact(j+1) * fact(n+1) * &
8781 & src_sk(j+1) / vscales(indj) / dst_si(n+1)
8784 tmp2 = ((-two)**(-l)) / &
8785 & fact(l+1)*fact(2*l+1)/ &
8786 & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
8787 & (work(j+n-l+1) / (z**l))
8791 res1 = res1 + tmp3*src_m(indj)
8801 tmp1 = vscales(indn+k) * ((-one)**(n)) * &
8802 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) * &
8803 & src_sk(j+1) / vscales(indj+k) / dst_si(n+1)
8806 tmp2 = ((-two)**(-l)) / &
8807 & fact(l+k+1)*fact(2*l+1)/ &
8808 & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
8809 & (work(j+n-l+1) / (z**l))
8813 res1 = res1 + tmp3*src_m(indj+k)
8814 res2 = res2 + tmp3*src_m(indj-k)
8816 dst_l(indn+k) = res1
8817 dst_l(indn-k) = res2
8831 tmp1 = vscales(indn) * ((-one)**(n)) * &
8832 & dble(2*j+1) * fact(j+1) * fact(n+1) * &
8833 & src_sk(j+1) / vscales(indj) / dst_si(n+1)
8836 tmp2 = ((-two)**(-l)) / &
8837 & fact(l+1)*fact(2*l+1)/ &
8838 & fact(l+1)/fact(l+1)/fact(j-l+1)/fact(n-l+1)* &
8839 & (work(j+n-l+1) / (z**l))
8843 res1 = res1 + tmp3*src_m(indj)
8845 dst_l(indn) = beta*dst_l(indn) + res1
8853 tmp1 = vscales(indn+k) * ((-one)**(n)) * &
8854 & dble(2*j+1) * fact(j-k+1) * fact(n+k+1) * &
8855 & src_sk(j+1) / vscales(indj+k) / dst_si(n+1)
8858 tmp2 = ((-two)**(-l)) / &
8859 & fact(l+k+1)*fact(2*l+1)/ &
8860 & fact(l+1)/fact(l-k+1)/fact(j-l+1)/fact(n-l+1)* &
8861 & (work(j+n-l+1) / (z**l))
8865 res1 = res1 + tmp3*src_m(indj+k)
8866 res2 = res2 + tmp3*src_m(indj-k)
8868 dst_l(indn+k) = beta*dst_l(indn+k) + res1
8869 dst_l(indn-k) = beta*dst_l(indn-k) + res2
8877 if (beta .eq. zero)
then
8882 do k = indj-j, indj+j
8883 dst_l(k) = tmp1 * src_m(k)
8893 do k = indj-j, indj+j
8894 dst_l(k) = beta*dst_l(k) + tmp1*src_m(k)
8929 & m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m)
8931 integer,
intent(in) :: pl, pm
8932 real(dp),
intent(in) :: c(3), src_r, dst_r, vscales((pm+pl+1)**2), &
8933 & m2l_ztranslate_adj_coef(pl+1, pl+1, pm+1), alpha, &
8934 & src_l((pl+1)*(pl+1)), beta
8936 real(dp),
intent(inout) :: dst_m((pm+1)*(pm+1))
8938 real(dp) :: work(6*max(pm, pl)**2 + 19*max(pm, pl) + 8)
8941 & m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m, work)
8973 & m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m, work)
8975 integer,
intent(in) :: pl, pm
8976 real(dp),
intent(in) :: c(3), src_r, dst_r, vscales((pm+pl+1)**2), &
8977 & m2l_ztranslate_adj_coef(pl+1, pl+1, pm+1), alpha, &
8978 & src_l((pl+1)*(pl+1)), beta
8980 real(dp),
intent(inout) :: dst_m((pm+1)*(pm+1))
8982 real(dp),
intent(out),
target :: &
8983 & work(6*max(pm, pl)**2 + 19*max(pm, pl) + 8)
8985 real(dp) :: rho, ctheta, stheta, cphi, sphi
8988 real(dp),
pointer :: tmp_ml(:), tmp_ml2(:), vcos(:), vsin(:)
8990 call carttosph(c, rho, ctheta, stheta, cphi, sphi)
8992 if (stheta .eq. 0)
then
8995 & m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m, work)
9002 tmp_ml(1:m) => work(n+1:n+m)
9004 tmp_ml2(1:m) => work(n+1:n+m)
9007 vcos => work(n+1:n+m)
9009 vsin => work(n+1:n+m)
9011 call trgev(cphi, sphi, p, vcos, vsin)
9019 & m2l_ztranslate_adj_coef, one, tmp_ml2, zero, tmp_ml, work)
Compile-time constants and definitions.
Harmonics-related core routines.
subroutine ylmscale(p, vscales, v4pi2lp1, vscales_rel)
Compute scaling factors of real normalized spherical harmonics.
subroutine fmm_m2p_adj(c, src_q, dst_r, p, vscales_rel, beta, dst_m)
Adjoint M2P operation.
subroutine fmm_l2p_adj_work(c, src_q, dst_r, p, vscales_rel, beta, dst_l, work)
Adjoint L2P.
subroutine fmm_m2m_ztranslate(z, src_r, dst_r, p, vscales, vcnk, alpha, src_m, beta, dst_m)
Direct M2M translation over OZ axis.
subroutine fmm_l2l_scale_adj(src_r, dst_r, p, alpha, src_l, beta, dst_l)
Adjoint scale L2L, when spherical harmonics are centered in the same point.
subroutine fmm_m2m_bessel_rotation_adj_work(c, src_sk, dst_sk, p, vscales, alpha, src_m, beta, dst_m, work, work_complex)
Direct M2M translation by 4 rotations and 1 translation.
subroutine fmm_m2l_bessel_rotation(c, src_r, dst_r, kappa, p, vscales, alpha, src_m, beta, dst_l)
Direct M2M translation by 4 rotations and 1 translation.
subroutine fmm_m2l_rotation(c, src_r, dst_r, pm, pl, vscales, m2l_ztranslate_coef, alpha, src_m, beta, dst_l)
Direct M2L translation by 4 rotations and 1 translation.
subroutine fmm_l2l_bessel_derivative_ztranslate_work(src_si, p, vscales, alpha, src_l, beta, dst_l)
Direct M2M translation over OZ axis.
subroutine fmm_m2m_ztranslate_work(z, src_r, dst_r, p, vscales, vcnk, alpha, src_m, beta, dst_m, work)
Direct M2M translation over OZ axis.
subroutine fmm_sph_rotate_oxz_work(p, ctheta, stheta, alpha, src, beta, dst, work)
Transform spherical harmonics in the OXZ plane.
subroutine fmm_m2m_ztranslate_adj(z, src_r, dst_r, p, vscales, vcnk, alpha, src_m, beta, dst_m)
Adjoint M2M translation over OZ axis.
subroutine fmm_m2m_bessel_derivative_ztranslate_work(src_sk, p, vscales, alpha, src_m, beta, dst_m)
Direct M2M translation over OZ axis.
subroutine fmm_l2l_rotation_adj(c, src_r, dst_r, p, vscales, vfact, alpha, src_l, beta, dst_l)
Adjoint L2L translation by 4 rotations and 1 translation.
subroutine fmm_m2l_ztranslate(z, src_r, dst_r, pm, pl, vscales, m2l_ztranslate_coef, alpha, src_m, beta, dst_l)
Direct M2L translation over OZ axis.
subroutine fmm_m2p_bessel_adj(c, src_q, dst_r, kappa, p, vscales, beta, dst_m)
Adjoint M2P operation.
subroutine fmm_m2m_rotation(c, src_r, dst_r, p, vscales, vcnk, alpha, src_m, beta, dst_m)
Direct M2M translation by 4 rotations and 1 translation.
subroutine fmm_l2l_scale(src_r, dst_r, p, alpha, src_l, beta, dst_l)
Scale L2L, when spherical harmonics are centered in the same point.
subroutine fmm_sph_rotate_oxz(p, ctheta, stheta, alpha, src, beta, dst)
Transform spherical harmonics in the OXZ plane.
subroutine fmm_m2l_ztranslate_adj_work(z, src_r, dst_r, pl, pm, vscales, m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m, work)
Adjoint M2L translation over OZ axis.
subroutine fmm_m2m_bessel_ztranslate(z, src_sk, dst_sk, p, vscales, alpha, src_m, beta, dst_m)
Direct M2M translation over OZ axis.
subroutine fmm_l2p(c, src_r, p, vscales, alpha, src_l, beta, dst_v)
Accumulate potential, induced by local spherical harmonics.
subroutine polleg_work(ctheta, stheta, p, vplm, work)
Compute associated Legendre polynomials.
subroutine fmm_l2l_ztranslate_work(z, src_r, dst_r, p, vscales, vfact, alpha, src_l, beta, dst_l, work)
Direct L2L translation over OZ axis.
subroutine fmm_l2p_adj(c, src_q, dst_r, p, vscales_rel, beta, dst_l)
Adjoint of L2P.
subroutine fmm_m2l_bessel_ztranslate(z, src_sk, dst_si, p, vscales, alpha, src_m, beta, dst_l)
Direct M2M translation over OZ axis.
subroutine fmm_p2m(c, src_q, dst_r, p, vscales, beta, dst_m)
Accumulate a multipole expansion induced by a particle of a given charge.
subroutine fmm_l2l_bessel_ztranslate_adj(z, src_si, dst_si, p, vscales, alpha, src_l, beta, dst_l)
Direct M2M translation over OZ axis.
subroutine fmm_m2m_scale(src_r, dst_r, p, alpha, src_m, beta, dst_m)
Scale M2M, when spherical harmonics are centered in the same point.
subroutine fmm_sph_rotate_oz_adj(p, vcos, vsin, alpha, src, beta, dst)
Rotate spherical harmonics around OZ axis in an opposite direction.
subroutine fmm_l2l_ztranslate_adj(z, src_r, dst_r, p, vscales, vfact, alpha, src_l, beta, dst_l)
Adjoint L2L translation over OZ axis.
subroutine fmm_l2l_bessel_ztranslate_work(z, src_si, dst_si, p, vscales, alpha, src_l, beta, dst_l, work, work_complex)
Direct M2M translation over OZ axis.
subroutine fmm_sph_rotate_oz_adj_work(p, vcos, vsin, alpha, src, beta, dst)
Rotate spherical harmonics around OZ axis in an opposite direction.
subroutine trgev(cphi, sphi, p, vcos, vsin)
Compute arrays of and .
subroutine fmm_m2p_adj_work(c, src_q, dst_r, p, vscales_rel, beta, dst_m, work)
Adjoint M2P operation.
subroutine carttosph(x, rho, ctheta, stheta, cphi, sphi)
Convert input cartesian coordinate into spherical coordinate.
subroutine fmm_l2l_bessel_rotation_adj_work(c, src_si, dst_si, p, vscales, alpha, src_l, beta, dst_l, work, work_complex)
Direct L2L translation by 4 rotations and 1 translation.
subroutine fmm_m2p_bessel_adj_work(c, src_q, dst_sk, p, vscales, beta, dst_m, work_complex, work)
Adjoint M2P operation.
subroutine fmm_m2m_bessel_rotation_work(c, src_sk, dst_sk, p, vscales, alpha, src_m, beta, dst_m, work, work_complex)
Direct M2M translation by 4 rotations and 1 translation.
subroutine ylmbas(x, rho, ctheta, stheta, cphi, sphi, p, vscales, vylm, vplm, vcos, vsin)
Compute all spherical harmonics up to a given degree at a given point.
subroutine fmm_l2l_rotation_work(c, src_r, dst_r, p, vscales, vfact, alpha, src_l, beta, dst_l, work)
Direct L2L translation by 4 rotations and 1 translation.
subroutine fmm_m2l_ztranslate_work(z, src_r, dst_r, pm, pl, vscales, m2l_ztranslate_coef, alpha, src_m, beta, dst_l, work)
Direct M2L translation over OZ axis.
subroutine polleg(ctheta, stheta, p, vplm)
Compute associated Legendre polynomials.
subroutine fmm_l2l_bessel_ztranslate_adj_work(z, src_si, dst_si, p, vscales, alpha, src_l, beta, dst_l, work, work_complex)
Direct M2M translation over OZ axis.
subroutine fmm_m2p_work(c, src_r, p, vscales_rel, alpha, src_m, beta, dst_v, work)
Accumulate potential, induced by multipole spherical harmonics.
subroutine fmm_l2l_rotation(c, src_r, dst_r, p, vscales, vfact, alpha, src_l, beta, dst_l)
Direct L2L translation by 4 rotations and 1 translation.
subroutine fmm_l2l_bessel_rotation(c, src_r, dst_r, kappa, p, vscales, alpha, src_l, beta, dst_l)
Direct L2L translation by 4 rotations and 1 translation.
subroutine fmm_m2l_rotation_adj(c, src_r, dst_r, pl, pm, vscales, m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m)
Adjoint M2L translation by 4 rotations and 1 translation.
subroutine fmm_sph_rotate_oz(p, vcos, vsin, alpha, src, beta, dst)
Rotate spherical harmonics around OZ axis.
subroutine fmm_m2m_rotation_work(c, src_r, dst_r, p, vscales, vcnk, alpha, src_m, beta, dst_m, work)
Direct M2M translation by 4 rotations and 1 translation.
subroutine fmm_sph_transform(p, r1, alpha, src, beta, dst)
Transform coefficients of spherical harmonics to a new cartesion system.
subroutine fmm_l2p_bessel(c, src_r, p, vscales, alpha, src_l, beta, dst_v)
Accumulate potential, induced by local spherical harmonics.
subroutine fmm_m2p_mat(c, r, p, vscales, mat)
Accumulate potentials, induced by each multipole spherical harmonic.
subroutine fmm_m2l_rotation_work(c, src_r, dst_r, pm, pl, vscales, m2l_ztranslate_coef, alpha, src_m, beta, dst_l, work)
Direct M2L translation by 4 rotations and 1 translation.
subroutine fmm_l2l_ztranslate_adj_work(z, src_r, dst_r, p, vscales, vfact, alpha, src_l, beta, dst_l, work)
Adjoint L2L translation over OZ axis.
subroutine fmm_l2l_bessel_rotation_adj(c, src_r, dst_r, kappa, p, vscales, alpha, src_l, beta, dst_l)
Direct L2L translation by 4 rotations and 1 translation.
subroutine fmm_l2l_rotation_adj_work(c, src_r, dst_r, p, vscales, vfact, alpha, src_l, beta, dst_l, work)
Adjoint L2L translation by 4 rotations and 1 translation.
subroutine fmm_m2m_bessel_ztranslate_adj(z, src_sk, dst_sk, p, vscales, alpha, src_m, beta, dst_m)
Direct M2M translation over OZ axis.
subroutine fmm_l2l_bessel_ztranslate(z, src_si, dst_si, p, vscales, alpha, src_l, beta, dst_l)
Direct M2M translation over OZ axis.
subroutine fmm_m2l_bessel_ztranslate_adj_work(z, src_sk, dst_si, p, vscales, alpha, src_m, beta, dst_l, work, work_complex)
Direct M2M translation over OZ axis.
subroutine fmm_m2m_bessel_ztranslate_adj_work(z, src_sk, dst_sk, p, vscales, alpha, src_m, beta, dst_m, work, work_complex)
Direct M2M translation over OZ axis.
subroutine fmm_m2m_ztranslate_adj_work(z, src_r, dst_r, p, vscales, vcnk, alpha, src_m, beta, dst_m, work)
Adjoint M2M translation over OZ axis.
subroutine fmm_l2l_ztranslate(z, src_r, dst_r, p, vscales, vfact, alpha, src_l, beta, dst_l)
Direct L2L translation over OZ axis.
subroutine fmm_sph_transform_work(p, r1, alpha, src, beta, dst, work)
Transform spherical harmonics to a new cartesion system of coordinates.
subroutine fmm_m2l_bessel_rotation_work(c, src_sk, dst_si, p, vscales, alpha, src_m, beta, dst_l, work, work_complex)
Direct M2M translation by 4 rotations and 1 translation.
subroutine fmm_m2l_bessel_rotation_adj(c, src_r, dst_r, kappa, p, vscales, alpha, src_m, beta, dst_l)
Direct M2M translation by 4 rotations and 1 translation.
subroutine fmm_l2l_bessel_rotation_work(c, src_si, dst_si, p, vscales, alpha, src_l, beta, dst_l, work, work_complex)
Direct L2L translation by 4 rotations and 1 translation.
subroutine fmm_m2l_rotation_adj_work(c, src_r, dst_r, pl, pm, vscales, m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m, work)
Adjoint M2L translation by 4 rotations and 1 translation.
subroutine fmm_m2p(c, src_r, p, vscales_rel, alpha, src_m, beta, dst_v)
Accumulate potential, induced by multipole spherical harmonics.
subroutine fmm_m2l_bessel_rotation_adj_work(c, src_sk, dst_si, p, vscales, alpha, src_m, beta, dst_l, work, work_complex)
Direct M2M translation by 4 rotations and 1 translation.
subroutine fmm_m2m_scale_adj(src_r, dst_r, p, alpha, src_m, beta, dst_m)
Adjoint scale M2M, when spherical harmonics are centered in the same point.
subroutine fmm_m2m_bessel_ztranslate_work(z, src_sk, dst_sk, p, vscales, alpha, src_m, beta, dst_m, work, work_complex)
Direct M2M translation over OZ axis.
subroutine fmm_m2l_bessel_ztranslate_work(z, src_sk, dst_si, p, vscales, alpha, src_m, beta, dst_l, work, work_complex)
Direct M2M translation over OZ axis.
subroutine fmm_m2l_bessel_ztranslate_adj(z, src_sk, dst_si, p, vscales, alpha, src_m, beta, dst_l)
Direct M2M translation over OZ axis.
subroutine fmm_m2p_bessel(c, src_r, p, vscales, alpha, src_m, beta, dst_v)
Accumulate potential, induced by multipole spherical harmonics.
subroutine fmm_constants(dmax, pm, pl, vcnk, m2l_ztranslate_coef, m2l_ztranslate_adj_coef)
Compute FMM-related constants.
subroutine fmm_p2m_work(c, src_q, dst_r, p, vscales, beta, dst_m, work)
Accumulate a multipole expansion induced by a particle of a given charge.
subroutine fmm_sph_rotate_oz_work(p, vcos, vsin, alpha, src, beta, dst)
Rotate spherical harmonics around OZ axis.
subroutine fmm_l2p_bessel_work(c, p, vscales, SI_ri, alpha, src_l, beta, dst_v, work_complex, work)
Accumulate Bessel potential, induced by local spherical harmonics.
subroutine fmm_m2m_bessel_rotation(c, src_r, dst_r, kappa, p, vscales, alpha, src_m, beta, dst_m)
Direct M2M translation by 4 rotations and 1 translation.
subroutine fmm_m2l_ztranslate_adj(z, src_r, dst_r, pl, pm, vscales, m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m)
Adjoint M2L translation over OZ axis.
subroutine fmm_m2m_rotation_adj_work(c, src_r, dst_r, p, vscales, vcnk, alpha, src_m, beta, dst_m, work)
Adjoint M2M translation by 4 rotations and 1 translation.
subroutine fmm_m2p_bessel_work(c, p, vscales, SK_ri, alpha, src_m, beta, dst_v, work_complex, work)
Accumulate Bessel potential, induced by multipole spherical harmonics.
subroutine fmm_m2m_bessel_rotation_adj(c, src_r, dst_r, kappa, p, vscales, alpha, src_m, beta, dst_m)
Direct M2M translation by 4 rotations and 1 translation.
subroutine fmm_l2p_work(c, src_r, p, vscales_rel, alpha, src_l, beta, dst_v, work)
Accumulate potential, induced by local spherical harmonics.
subroutine fmm_m2m_rotation_adj(c, src_r, dst_r, p, vscales, vcnk, alpha, src_m, beta, dst_m)
Adjoint M2M translation by 4 rotations and 1 translation.