1c program test 2c implicit none 3c 4c integer nion,n,na,nqm,nb,lr 5c real*8 a_in,b_in,aqm,da,a_out,b_out 6c 7c write(*,*) "Enter a,b,aqm,nion,lr:" 8c read(*,*) a_in,b_in,aqm,nion,lr 9c 10c call FMM_Generate_Box(a_in,b_in,aqm,nion,lr, 11c > n,na,nqm,nb,da,a_out,b_out) 12c 13c write(*,*) 14c write(*,*) "da,n,n*da = ",n,da,da*n 15c write(*,*) "a_in,a_out,na = ",na,a_in,a_out 16c write(*,*) "b_in,b_out,nb = ",nb,b_in,b_out 17c write(*,*) "aqm,nqm = ",nqm,aqm 18c write(*,*) 19c stop 20c end 21 22* ********************************************** 23* * * 24* * FMM_start * 25* * * 26* ********************************************** 27* 28 subroutine FMM_start() 29 implicit none 30 31#include "util.fh" 32#include "stdio.fh" 33#include "fmm.fh" 34 35* **** local variables **** 36 integer MASTER,taskid,np 37 parameter (MASTER=0) 38 39 logical oprint 40 integer l,m,i,j,k,lm 41 42* **** external functions **** 43 logical control_fmm,control_print 44 integer control_fmm_lmax,control_fmm_lr,ion_nion,control_version 45 external control_fmm,control_print 46 external control_fmm_lmax,control_fmm_lr,ion_nion,control_version 47 48 fmm = control_fmm().and.(control_version().eq.4) 49 if (fmm) then 50 lmax = control_fmm_lmax() 51 lr = control_fmm_lr() 52 levels = nint(dlog(dble(ion_nion()))/dlog(8.0d0)) - lr 53 nboxes = 8**levels 54 tboxes = 0 55 do l=0,levels 56 tboxes = tboxes + 8**l 57 end do 58 59 call Parallel_taskid(taskid) 60 call Parallel_np(np) 61 oprint = (taskid.eq.MASTER).and.control_print(print_medium) 62 63c value = BA_alloc_get(mt_dbl,26*(lmax+1)**2,'first_tlm_fmm', 64c > first_tlm_fmm(2),first_tlm_fmm(1)) 65c value = value.and. 66c > BA_alloc_get(mt_dbl,98*(lmax+1)**2,'second_tlm_fmm', 67c > second_tlm_fmm(2),second_tlm_fmm(1)) 68c value = value.and. 69c > BA_alloc_get(mt_dbl,7216,'interaction_tlm_fmm', 70c > interaction_tlm_fmm(2),interaction_tlm_fmm(1)) 71c 72c* **** first neigbors of the current box **** 73c* **** first_count = 3*3*3 - 1*1*1 = 27-1 = 26 **** 74c first_count = 0 75c do k=-1,1 76c do j=-1,1 77c do i=-1,1 78c if ((i.ne.0).or.(j.ne.0).or.(k.eq.0)) then 79c do l=0,fmm_lmax 80c do m=-l,l 81c dbl_mb(first_tlm_fmm(1)+first_count) 82c > = Tesseral3_lm(l,m,dble(i),dble(j),dble(k)) 83c first_count = first_count + 1 84c end do 85c end do 86c end if 87c end do 88c end do 89c end do 90c 91c* **** second nearest neigbors of the current box **** 92c* **** second_count = 5*5*5 - 3*3*3 = 125-27 = 98 **** 93c second_count = 0 94c do k=-2,2 95c do j=-2,2 96c do i=-2,2 97c if ((i*i+j*j+k*k).gt.2) then 98c do l=0,fmm_lmax 99c do m=-l,l 100c dbl_mb(second_tlm_fmm(1)+second_count) 101c > = Tesseral3_lm(l,m,dble(i),dble(j),dble(k)) 102c second_count = second_count + 1 103c end do 104c end do 105c end if 106c end do 107c end do 108c end do 109 110c* **** interaction neigbors - set of boxes which are children of **** 111c* **** the nearest and second nearest neighbors of the current boxes **** 112c* **** parent which are not nearest or second nearest neighbors of **** 113c* **** the current box **** 114c* **** interaction_count=8*(8*5*5*5-98)=8*(8*125-98)=8*(1000-98)=8*902=7216 **** 115 116 if (oprint) then 117 write(luout,100) lmax,lr, 118 > levels, 119 > nboxes,tboxes 120 end if 121 end if 122 return 123 100 format(/" FMM - Fast Multipole Algorithm", 124 > /" - fmm_lmax = ",i3," fmm_lr = ",i3, 125 > /" - number of refinement levels = ",i2, 126 > /" - nboxes (finest level) = ",i8, 127 > /" - nboxes (total) = ",i8/) 128 end 129 130 131* ********************************************** 132* * * 133* * FMM_end * 134* * * 135* ********************************************** 136 137 subroutine FMM_end() 138 implicit none 139 140#include "fmm.fh" 141 142 if (fmm) then 143 end if 144 return 145 end 146 147 148 logical function FMM_fmm() 149 implicit none 150 151#include "fmm.fh" 152 153 FMM_fmm = fmm 154 return 155 end 156 157 integer function FMM_lmax() 158 implicit none 159 160#include "fmm.fh" 161 162 FMM_lmax = lmax 163 return 164 end 165 166* ********************************************** 167* * * 168* * FMM_Setup_Tree * 169* * * 170* ********************************************** 171 172* This routine generates a box gridding and tree for FMM calculations. 173* 174* Entry - nion - number of ions 175* rion - atom positions 176* levels - number of tree levels 177* 178* Exit - a_out - box length 179* ion_boxindx - box indexes for the ions 180 181 subroutine FMM_Setup_Tree(nion,rion,levels,a_out,ion_boxindx) 182 implicit none 183 integer nion 184 real*8 rion(3,*) 185 integer levels 186 187 real*8 a_out 188 integer ion_boxindx(3,*) 189 190* **** local variables **** 191 integer ii,n,n2 192 real*8 aover2 193 194 aover2 = 0.0d0 195 do ii=1,nion 196 if (dabs(rion(1,ii)).gt.aover2) aover2=dabs(rion(1,ii)) 197 if (dabs(rion(2,ii)).gt.aover2) aover2=dabs(rion(2,ii)) 198 if (dabs(rion(3,ii)).gt.aover2) aover2=dabs(rion(3,ii)) 199 end do 200 aover2 = aover2 + 0.5d0 201 a_out = 2*aover2 202 203 n = 2**levels 204 n2 = n*n 205 do ii=1,nion 206 ion_boxindx(1,ii) = idint(n*(rion(1,ii)+aover2)/a_out) 207 ion_boxindx(2,ii) = idint(n*(rion(2,ii)+aover2)/a_out) 208 ion_boxindx(3,ii) = idint(n*(rion(3,ii)+aover2)/a_out) 209 end do 210 211 return 212 end 213 214* ********************************************** 215* * * 216* * FMM_Setup_Mlm * 217* * * 218* ********************************************** 219 220 subroutine FMM_Setup_Mlm(nion,katm,rion,qatm, 221 > levels,a,ion_boxindx, 222 > lmax,Mlm) 223 implicit none 224 integer nion,katm(*) 225 real*8 rion(3,*) 226 real*8 qatm(*) 227 integer levels 228 real*8 a 229 integer ion_boxindx(3,*) 230 integer lmax 231 complex*16 Mlm((lmax+1)**2,*) 232 233* **** local variables **** 234 integer taskid,np 235 integer ii,i,j,k,n,n2,ibox,idx,l,m 236 real*8 aover2,da,xc,yc,zc,q 237 238* **** external functions **** 239 complex*16 rSphericalHarmonic3 240 external rSphericalHarmonic3 241 242 call Parallel_np(np) 243 call Parallel_taskid(taskid) 244 245 n = 2**levels 246 n2 = n*n 247 da = a/dble(a) 248 aover2 = a/2.0d0 249 call dcopy(2*((lmax+1)**2)*(n*n*n),0.0d0,0,Mlm,1) 250 do ii=1,nion 251 if (mod(ii,np).eq.taskid) then 252 q = qatm(katm(ii)) 253 xc = -aover2 + (ion_boxindx(1,ii)+0.5d0)*da 254 yc = -aover2 + (ion_boxindx(2,ii)+0.5d0)*da 255 zc = -aover2 + (ion_boxindx(3,ii)+0.5d0)*da 256 ibox = 1 + ion_boxindx(1,ii) 257 > + ion_boxindx(2,ii)*n 258 > + ion_boxindx(3,ii)*n2 259 idx = 1 260 do l=0,lmax 261 do m=-l,l 262 Mlm(idx,ibox) = Mlm(idx,ibox) 263 > + q*rSphericalHarmonic3(l,-m, 264 > xc-rion(1,ii), 265 > yc-rion(2,ii), 266 > zc-rion(3,ii)) 267 idx = idx + 1 268 end do 269 end do 270 end if 271 end do 272 call Parallel_Vector_SumAll(((lmax+1)**2)*(n*n*n),Mlm) 273 274 return 275 end 276 277* ********************************************** 278* * * 279* * FMM_Generate_Box * 280* * * 281* ********************************************** 282* 283* This routine generates a box gridding for FMM calculations. 284* 285* Entry - a_in - guess for start of box 286* b_in - guess for end of box 287* aqm - simple cubic side length 288* nion - number of ions 289* lr - level reduction factor 290* 291* Exit - levels - number of tree levels 292* n - number of boxes from a to b 293* na - number of boxes from a to (-aqm/2) 294* nqm - number of boxes from (-aqm/2) to (aqm/2) 295* nb - number of boxes from (aqm/2) to b 296* da - length of small boxes 297* a_out - box start 298* b_out - box end 299* 300 subroutine FMM_Generate_Box(a_in,b_in,aqm,nion,lr, 301 > levels,n,na,nqm,nb,da,a_out,b_out) 302 implicit none 303 real*8 a_in,b_in,aqm 304 integer nion,lr 305 integer levels,n,na,nqm,nb 306 real*8 da,a_out,b_out 307 308 levels = idint(dlog(1.0d0*nion)/dlog(8.0d0))-lr 309 n = 2**levels 310 311 nqm = idint(n*aqm/(b_in-a_in)) 312 da = aqm/dble(nqm) 313 314 na = idint(n*(-0.5d0*aqm-a_in)/(b_in-a_in)) 315 nb = n - nqm - na 316 a_out = -0.5d0*aqm - na*da 317 b_out = 0.5d0*aqm + nb*da 318 return 319 end 320 321 322 323* ********************************************** 324* * * 325* * FMM_rion_Llm * 326* * * 327* ********************************************** 328* 329 subroutine FMM_rion_Llm(lmax,nion,nion_qm,katm,rion,qatm, 330 > lmbda0,rsphere2,Llm) 331 implicit none 332 integer lmax 333 integer nion,nion_qm,katm(*) 334 real*8 rion(3,*),qatm(*),lmbda0,rsphere2,Llm(*) 335 336* **** local variables **** 337 integer taskid,np 338 integer ii,l,m,idx 339 real*8 fourpi,r2,lmbda 340 341* **** external functions **** 342 real*8 Tesseral3_lm_over_rl 343 external Tesseral3_lm_over_rl 344 345 call Parallel_np(np) 346 call Parallel_taskid(taskid) 347 348 fourpi = 16.0d0*datan(1.0d0) 349 !call dcopy((lmax+1)**2,0.0d0,0,Llm,1) 350 call Parallel_shared_vector_zero(.true.,(lmax+1)**2,Llm) 351 do ii=1,nion 352 if (mod(ii,np).eq.taskid) then 353 r2 = rion(1,ii)**2+rion(2,ii)**2+rion(3,ii)**2 354 if (r2.gt.rsphere2) then 355 if (ii.gt.nion_qm) then 356 lmbda = lmbda0 357 else 358 lmbda = 1.0d0 359 end if 360 !idx = 1 361!$OMP DO 362 do l=0,lmax 363 do m=-l,l 364 idx = l*l+l+m+1 365 Llm(idx) = Llm(idx) 366 > + (fourpi/dble(2*l+1)) 367 > *qatm(katm(ii))*lmbda 368 > *Tesseral3_lm_over_rl(l,m,rion(1,ii), 369 > rion(2,ii), 370 > rion(3,ii)) 371 !idx = idx + 1 372 end do 373 end do 374!$OMP END DO 375 end if 376 end if 377 end do 378 call Parallel_Vector_SumAll((lmax+1)**2,Llm) 379 return 380 end 381 382* ********************************************** 383* * * 384* * FMM_rion_Llm_mm * 385* * * 386* ********************************************** 387* 388 subroutine FMM_rion_Llm_mm(lmax,nion,nion_qm,katm,rion,qatm, 389 > rsphere2,Llm) 390 implicit none 391 integer lmax 392 integer nion,nion_qm,katm(*) 393 real*8 rion(3,*),qatm(*),lmbda0,rsphere2,Llm(*) 394 395* **** local variables **** 396 integer taskid,np 397 integer ii,l,m,idx 398 real*8 fourpi,r2 399 400* **** external functions **** 401 real*8 Tesseral3_lm_over_rl 402 external Tesseral3_lm_over_rl 403 404 call Parallel_np(np) 405 call Parallel_taskid(taskid) 406 407 fourpi = 16.0d0*datan(1.0d0) 408 !call dcopy((lmax+1)**2,0.0d0,0,Llm,1) 409 call Parallel_shared_vector_zero(.true.,(lmax+1)**2,Llm) 410 do ii=nion_qm+1,nion 411 if (mod(ii,np).eq.taskid) then 412 r2 = rion(1,ii)**2+rion(2,ii)**2+rion(3,ii)**2 413 if (r2.gt.rsphere2) then 414 415 !idx = 1 416!$OMP DO 417 do l=0,lmax 418 do m=-l,l 419 idx = l*l+l+m+1 420 Llm(idx) = Llm(idx) 421 > + (fourpi/dble(2*l+1)) 422 > *qatm(katm(ii)) 423 > *Tesseral3_lm_over_rl(l,m,rion(1,ii), 424 > rion(2,ii), 425 > rion(3,ii)) 426 !idx = idx + 1 427 end do 428 end do 429!$OMP END DO 430 end if 431 end if 432 end do 433 call Parallel_Vector_SumAll((lmax+1)**2,Llm) 434 return 435 end 436 437 438 439 440* ********************************************** 441* * * 442* * FMM_fion_Mlm * 443* * * 444* ********************************************** 445* 446 subroutine FMM_fion_Mlm(lmax,nion,nion_qm,katm,rion,qatm, 447 > lmbda0,rsphere2, 448 > Mlm,fion) 449 implicit none 450 integer lmax 451 integer nion,nion_qm,katm(*) 452 real*8 rion(3,*),qatm(*),lmbda0,rsphere2,Mlm(*),fion(3,nion) 453 454* **** local variables **** 455 integer taskid,np 456 integer ii,l,m,idx 457 real*8 fourpi,r2,coef,dtx,dty,dtz,lmbda 458 459 460 call Parallel_np(np) 461 call Parallel_taskid(taskid) 462 463 fourpi = 16.0d0*datan(1.0d0) 464 do ii=1,nion 465 if (mod(ii,np).eq.taskid) then 466 r2 = rion(1,ii)**2+rion(2,ii)**2+rion(3,ii)**2 467 if (r2.gt.rsphere2) then 468 if (ii.gt.nion_qm) then 469 lmbda = lmbda0 470 else 471 lmbda = 1.0d0 472 end if 473c idx = 1 474!$OMP DO REDUCTION(+:fion) 475 do l=0,lmax 476 coef = (fourpi/dble(2*l+1))*qatm(katm(ii)) 477 do m=-l,l 478 idx = l*l+l+m+1 479 call dTesseral3_lm_over_rl(l,m,rion(1,ii), 480 > rion(2,ii), 481 > rion(3,ii), 482 > dtx,dty,dtz) 483 fion(1,ii) = fion(1,ii) + dtx*coef*Mlm(idx)*lmbda 484 fion(2,ii) = fion(2,ii) + dty*coef*Mlm(idx)*lmbda 485 fion(3,ii) = fion(3,ii) + dtz*coef*Mlm(idx)*lmbda 486c idx = idx + 1 487 end do 488 end do 489!$OMP END DO 490 end if 491 end if 492 end do 493 return 494 end 495 496 497 498c* ********************************************** 499c* * * 500c* * FMM_rgrid_Mlm * 501c* * * 502c* ********************************************** 503c 504c subroutine FMM_rgrid_Mlm(lmax,ngrid,rgrid,Mlm) 505c implicit none 506c integer lmax 507c integer ngrid 508c real*8 rgrid(3,*) 509c real*8 Mlm(*) 510c 511c* **** local variables **** 512c integer i,l,m,idx 513c real*8 x,y,z 514c 515c* **** external functions **** 516c real*8 Tesseral3_lm_rl 517c external Tesseral3_lm_rl 518c 519c do i=1,ngrid 520c x = rgrid(1,i) 521c y = rgrid(2,i) 522c z = rgrid(3,i) 523c vl(i) = 0.0d0 524c idx = 1 525c do l=0,lmax 526c do m=-l,l 527c vl(i) = vl(i) + Llm(idx)*Tesseral3_lm_over_rl(l,m,x,y,z) 528c idx = idx + 1 529c end do 530c end do 531c end do 532c return 533c end 534 535* *************************************** 536* * * 537* * FMM_A_over_i * 538* * * 539* *************************************** 540 complex*16 function FMM_A_over_i(l,m) 541 implicit none 542 integer l,m 543 544* **** local variables **** 545 integer i,mod_m 546 real*8 A 547 complex*16 img 548 549 550 img = dcmplx(0.0d0,1.0d0) 551 mod_m = abs(m) 552 A = 1 553 do i=1,(l-mod_m) 554 A = A*dble(i) 555 end do 556 A = A*A 557 do i=(l-mod_m+1),(l+mod_m) 558 A = A*dble(i) 559 end do 560 A = ((-1)**l)/dsqrt(A) 561 562 FMM_A_over_i = A/img**mod_m 563 return 564 end 565 566* ********************************************** 567* * * 568* * FMM_Translate_Mlm * 569* * * 570* ********************************************** 571 572 subroutine FMM_Translate_Mlm(lmax,Olm,Nlm,x,y,z) 573 implicit none 574 integer lmax 575 complex*16 Olm(*) 576 complex*16 Nlm(*) 577 real*8 x,y,z 578 579#include "bafdecls.fh" 580#include "errquit.fh" 581 582 583* **** local variables **** 584 logical value 585 integer l,m,idx,n 586 integer rYlm(2),A(2),B(2),AB(2),tmp22(2),tmp42(2) 587 588* **** external functions **** 589 complex*16 rSphericalHarmonic3 590 external rSphericalHarmonic3 591 592 n = (2*lmax+2)*(4*lmax+2) 593 594* **** allocate stack **** 595 value = BA_push_get(mt_dcpl,(lmax+1)**2,'rYlm',rYlm(2),rYlm(1)) 596 value=value.and.BA_push_get(mt_dcpl,n,'A', A(2), A(1)) 597 value=value.and.BA_push_get(mt_dcpl,n,'B', B(2), B(1)) 598 value=value.and.BA_push_get(mt_dcpl,n,'AB',AB(2),AB(1)) 599 value=value.and. 600 > BA_push_get(mt_dcpl,(2*lmax+2),'tmp22',tmp22(2),tmp22(1)) 601 value=value.and. 602 > BA_push_get(mt_dcpl,(4*lmax+2),'tmp42',tmp42(2),tmp42(1)) 603 if (.not.value) 604 > call errquit('FMM_Translate_Mlm:out of stack',0,MA_ERR) 605 606 call dcffti(2*lmax+2,dcpl_mb(tmp22(1))) 607 call dcffti(4*lmax+2,dcpl_mb(tmp42(1))) 608 609 idx = 0 610 do l=0,lmax 611 do m=-l,l 612 dcpl_mb(rYlm(1)+idx) = rSphericalHarmonic3(l,-m,x,y,z) 613 end do 614 end do 615 616 call FMM_Translate_Mlm_sub1(lmax,Olm,dcpl_mb(rYlm(1)), 617 > dcpl_mb(A(1)),dcpl_mb(B(1))) 618 619 call FMM_Translate_sub2(lmax,dcpl_mb(A(1)),dcpl_mb(B(1)), 620 > dcpl_mb(AB(1)), 621 > dcpl_mb(tmp22(1)),dcpl_mb(tmp42(1))) 622 623 call FMM_Translate_Mlm_sub3(lmax,dcpl_mb(AB(1)),Nlm) 624 625 626* **** deallocate stack **** 627 value = BA_pop_stack(tmp42(2)) 628 value = value.and.BA_pop_stack(tmp22(2)) 629 value = value.and.BA_pop_stack(AB(2)) 630 value = value.and.BA_pop_stack(B(2)) 631 value = value.and.BA_pop_stack(A(2)) 632 value = value.and.BA_pop_stack(rYlm(2)) 633 if (.not.value) 634 > call errquit('FMM_Translate_Mlm:popping stack',0,MA_ERR) 635 return 636 end 637 638* *************************************** 639* * * 640* * FMM_Translate_Mlm_sub1 * 641* * * 642* *************************************** 643 subroutine FMM_Translate_Mlm_sub1(lmax,Olm,rYlm,A,B) 644 implicit none 645 integer lmax 646 complex*16 Olm(*),rYlm(*) 647 complex*16 A(2*lmax+2,4*lmax+2) 648 complex*16 B(2*lmax+2,4*lmax+2) 649 650 651* **** local variables **** 652 integer idx,i,j,l,m 653 complex*16 ctmp 654 655* **** external functions **** 656 complex*16 FMM_A_over_i 657 external FMM_A_over_i 658 659 call dcopy(2*(2*lmax+2)*(4*lmax+2),0.0d0,0,A,1) 660 call dcopy(2*(2*lmax+2)*(4*lmax+2),0.0d0,0,B,1) 661 idx = 1 662 j = 1 663 do l=0,lmax 664 i=1 665 do m=-l,l 666 ctmp = FMM_A_over_i(l,m) 667 A(i,j) = Olm(idx) *ctmp 668 B(i,j) = rYlm(idx)*ctmp 669 670 i = i + 1 671 idx = idx + 1 672 end do 673 j = j + 1 674 end do 675 return 676 end 677 678* *************************************** 679* * * 680* * FMM_Translate_sub2 * 681* * * 682* *************************************** 683 subroutine FMM_Translate_sub2(lmax,A,B,AB,tmp22,tmp42) 684 implicit none 685 integer lmax 686 complex*16 A(2*lmax+2,4*lmax+2) 687 complex*16 B(2*lmax+2,4*lmax+2) 688 complex*16 AB(2*lmax+2,4*lmax+2) 689 complex*16 tmp22(*),tmp42(*) 690 691* **** local variables **** 692 integer i,j,n 693 694 n = (2*lmax+2)*(4*lmax+2) 695 do j=1,(4*lmax+2) 696 call dcfftf(2*lmax+2,A(1,j),tmp22) 697 call dcfftf(2*lmax+2,B(1,j),tmp22) 698 end do 699 do i=1,(2*lmax+2) 700 call dcfftf(4*lmax+2,A(i,1),tmp42) 701 call dcfftf(4*lmax+2,B(i,1),tmp42) 702 end do 703 do j=1,(4*lmax+2) 704 do i=1,(2*lmax+2) 705 AB(i,j) = A(i,j)*B(i,j) 706 end do 707 end do 708 709 do i=1,(2*lmax+2) 710 call dcfftb(4*lmax+2,AB(i,1),tmp42) 711 end do 712 do j=1,(4*lmax+2) 713 call dcfftb(2*lmax+2,AB(1,j),tmp22) 714 end do 715 call dscal(2*n,1.0d0/dble(n),AB,1) 716 717 return 718 end 719 720* *************************************** 721* * * 722* * FMM_Translate_Mlm_sub3 * 723* * * 724* *************************************** 725 subroutine FMM_Translate_Mlm_sub3(lmax,AB,Mlm) 726 implicit none 727 integer lmax 728 complex*16 AB(2*lmax+2,4*lmax+2) 729 complex*16 Mlm(*) 730 731* **** local variables **** 732 integer idx,i,j,l,m 733 734* **** external functions **** 735 complex*16 FMM_A_over_i 736 external FMM_A_over_i 737 738* Jeff added these because they are implicit with Intel Fortran 739* and either this code is never used, is fine to give undefined behavior 740* or assumes initialization to zero. 741 idx = 0 742 j = 0 743 i = 0 744 745 do l=0,lmax 746 do m=-l,l 747 Mlm(idx) = AB(i,j)/FMM_A_over_i(l,m) 748 i = i + 1 749 idx = idx + 1 750 end do 751 j = j + 1 752 end do 753 return 754 end 755 756 757* ********************************************** 758* * * 759* * FMM_Translate_Llm * 760* * * 761* ********************************************** 762 763 subroutine FMM_Translate_Llm(lmax,Olm,Nlm,x,y,z) 764 implicit none 765 integer lmax 766 complex*16 Olm(*) 767 complex*16 Nlm(*) 768 real*8 x,y,z 769 770#include "bafdecls.fh" 771#include "errquit.fh" 772 773 774* **** local variables **** 775 logical value 776 integer l,m,idx,n 777 integer rYlm(2),A(2),B(2),AB(2),tmp22(2),tmp42(2) 778 779* **** external functions **** 780 complex*16 rSphericalHarmonic3 781 external rSphericalHarmonic3 782 783 n = (2*lmax+2)*(4*lmax+2) 784 785* **** allocate stack **** 786 value = BA_push_get(mt_dcpl,(lmax+1)**2,'rYlm',rYlm(2),rYlm(1)) 787 value=value.and.BA_push_get(mt_dcpl,n,'A', A(2), A(1)) 788 value=value.and.BA_push_get(mt_dcpl,n,'B', B(2), B(1)) 789 value=value.and.BA_push_get(mt_dcpl,n,'AB',AB(2),AB(1)) 790 value=value.and. 791 > BA_push_get(mt_dcpl,2*lmax+2,'tmp22',tmp22(2),tmp22(1)) 792 value=value.and. 793 > BA_push_get(mt_dcpl,4*lmax+2,'tmp42',tmp42(2),tmp42(1)) 794 if (.not.value) 795 > call errquit('FMM_Translate_Llm:out of stack',0,MA_ERR) 796 797 call dcffti(2*lmax+2,dcpl_mb(tmp22(1))) 798 call dcffti(4*lmax+2,dcpl_mb(tmp42(1))) 799 800 idx = 0 801 do l=0,lmax 802 do m=-l,l 803 dcpl_mb(rYlm(1)+idx) = rSphericalHarmonic3(l,m,x,y,z) 804 end do 805 end do 806 807 call FMM_Translate_Llm_sub1(lmax,Olm,dcpl_mb(rYlm(1)), 808 > dcpl_mb(A(1)),dcpl_mb(B(1))) 809 810 call FMM_Translate_sub2(lmax,dcpl_mb(A(1)),dcpl_mb(B(1)), 811 > dcpl_mb(AB(1)), 812 > dcpl_mb(tmp22(1)),dcpl_mb(tmp42(1))) 813 814 call FMM_Translate_Llm_sub3(lmax,dcpl_mb(AB(1)),Nlm) 815 816 817* **** deallocate stack **** 818 value = BA_pop_stack(tmp42(2)) 819 value = value.and.BA_pop_stack(tmp22(2)) 820 value = value.and.BA_pop_stack(AB(2)) 821 value = value.and.BA_pop_stack(B(2)) 822 value = value.and.BA_pop_stack(A(2)) 823 value = value.and.BA_pop_stack(rYlm(2)) 824 if (.not.value) 825 > call errquit('FMM_Translate_Llm:popping stack',0,MA_ERR) 826 return 827 end 828 829 830* *************************************** 831* * * 832* * FMM_Translate_Llm_sub1 * 833* * * 834* *************************************** 835 subroutine FMM_Translate_Llm_sub1(lmax,Olm,rYlm,A,B) 836 implicit none 837 integer lmax 838 complex*16 Olm(*),rYlm(*) 839 complex*16 A(2*lmax+2,4*lmax+2) 840 complex*16 B(2*lmax+2,4*lmax+2) 841 842 843* **** local variables **** 844 integer idx,i,j,l,m,sgn 845 complex*16 ctmp 846 847* **** external functions **** 848 complex*16 FMM_A_over_i 849 external FMM_A_over_i 850 851 call dcopy(2*(2*lmax+2)*(4*lmax+2),0.0d0,0,A,1) 852 call dcopy(2*(2*lmax+2)*(4*lmax+2),0.0d0,0,B,1) 853 idx = 1 854 j = 1 855 sgn = 1 856 do l=0,lmax 857 i=1 858 do m=-l,l 859 ctmp = FMM_A_over_i(l,m) 860 A(i,j) = sgn*Olm(idx)/ctmp 861 B(i,j) = rYlm(idx)*ctmp 862 863 i = i + 1 864 idx = idx + 1 865 end do 866 j = j + 1 867 sgn = -sgn 868 end do 869 return 870 end 871 872 873 874* *************************************** 875* * * 876* * FMM_Translate_Llm_sub3 * 877* * * 878* *************************************** 879 subroutine FMM_Translate_Llm_sub3(lmax,AB,Llm) 880 implicit none 881 integer lmax 882 complex*16 AB(2*lmax+2,4*lmax+2) 883 complex*16 Llm(*) 884 885* **** local variables **** 886 integer idx,i,j,l,m,sgn 887 888* **** external functions **** 889 complex*16 FMM_A_over_i 890 external FMM_A_over_i 891 892* Jeff: idx was not initialized before and defaults to zero with Intel, 893* which appears to be the primary compiler for testing. 894 idx = 0 895 896 sgn = 1 897 j=1 898 do l=0,lmax 899 i=1 900 do m=-l,l 901 Llm(idx) = sgn*AB(i,j)/FMM_A_over_i(l,m) 902 i = i + 1 903 idx = idx + 1 904 end do 905 j = j + 1 906 sgn = -sgn 907 end do 908 return 909 end 910 911 912* *************************************** 913* * * 914* * FMM_reduce_levels * 915* * * 916* *************************************** 917* 918* This function translates index values between different levels of an octree grid 919 920 integer function FMM_reduce_levels(nl,l,lr) 921 implicit none 922 integer nl,l,lr 923 924 FMM_reduce_levels = ishft(ishft(ibits(nl,2*l,l),-lr),2*(l-lr)) 925 > + ishft(ishft(ibits(nl,l,l),-lr),(l-lr)) 926 > + ishft(ibits(nl,0,l),-lr) 927 return 928 end 929 930 931