1! 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt. 6! See Docs/Contributors.txt for a list of contributors. 7! 8module m_delk 9 10 implicit none 11 12 private 13 14 public :: delk 15 16contains 17 18 subroutine delk( iexpikr, no, np, dvol, nvmax, & 19 numVs, listVsptr, listVs, & 20 nuo, nuotot, iaorb, iphorb, isa ) 21 22! ******************************************************************** 23! Finds the matrix elements of a plane wave 24! < \phi_{\mu} | exp^( iexpikr * i * \vec{k} \cdot \vec{r} ) | \phi_{\nu} > 25! First version written by J. Junquera in Feb. 2008 26! Adapted from an existing version of vmat after the parallelization 27! designed by BSC in July 2011. 28! *********************** INPUT ************************************** 29! integer iexpikr : Prefactor of the dot product between the 30! the k-vector and the position-vector in exponent. 31! it might be +1 or -1 32! integer no : Number of basis orbitals 33! integer np : Number of columns in C (local) 34! real*8 dvol : Volume per mesh point 35! integer nvmax : First dimension of listV and Vs, and maxim. 36! number of nonzero elements in any row of delkmat 37! integer numVs(nuo) : Number of non-zero elements in a row of delkmat 38! integer listVsptr(nuo) : Pointer to the start of rows in listVs 39! integer listVs(nvmax) : List of non-zero elements of delkmat 40! integer iaorb(*) : Pointer to atom to which orbital belongs 41! integer iphorb(*) : Orbital index within each atom 42! integer isa(*) : Species index of all atoms 43! ***************************** OUTPUT ******************************* 44! complex*16 delkmat(nvmax) : Value of nonzero elements in each row 45! of the matrix elements of exp(i*\vec{k}\vec{r}) 46! this variable is defined in the module 47! m_dimensionsefield (file m_dimefield.f90) 48! ********************************************************************* 49 50! Modules 51 use precision, only: dp, grid_p 52 use atmfuncs, only: rcut, all_phi 53 use atm_types, only: nsmax=>nspecies 54 use atomlist, only: indxuo, indxua 55 use siesta_geom, only: xa 56 use listsc_module, only: listsc 57 use mesh, only: dxa, nsp, xdop, xdsp, ne, nem 58 use mesh, only: cmesh, ipa, idop, nmsc, iatfold 59 use mesh, only: meshLim 60 use meshdscf, only: matrixMtoOC 61 use meshdscf, only: needdscfl, listdl, numdl, nrowsdscfl, listdlptr 62 use meshphi, only: directphi, endpht, lstpht, listp2, phi 63 use parallel, only: Nodes, node 64 use alloc, only: re_alloc, de_alloc 65 use parallelsubs, only: GlobalToLocalOrb 66 use m_planewavematrixvar, only: delkmat, wavevector 67#ifdef MPI 68 use mpi_siesta 69#endif 70#ifdef _OPENMP 71 use omp_lib 72#endif 73 74! Argument types and dimensions 75 integer :: iexpikr, no, np, nvmax, nuo, nuotot 76 integer :: iaorb(*), iphorb(*), isa(*), numVs(nuo) 77 integer :: listVsptr(nuo), listVs(nvmax) 78 real(dp) :: dvol 79! Internal variables and arrays 80 integer, parameter :: minloc = 1000 ! Min buffer size 81 integer, parameter :: maxoa = 100 ! Max # of orb/atom 82 integer :: i, ia, ic, ii, ijl, il, imp, ind, iop 83 integer :: ip, iphi, io, is, isp, iu, iul 84 integer :: ix, j, jc, jl, last, lasta, lastop 85 integer :: maxloc, maxloc2, nc, nlocal, nphiloc 86 integer :: nvmaxl, triang, lenx, leny, lenz,lenxy 87 logical :: ParallelLocal 88 real(dp) :: Vij(2), r2sp, dxsp(3), VClocal(2,nsp) 89 90 integer, pointer :: ilc(:), ilocal(:), iorb(:) 91 real(dp), pointer :: DscfL(:,:), t_DscfL(:,:,:), Clocal(:,:) 92 real(dp), pointer :: Vlocal(:,:), phia(:,:), r2cut(:) 93 integer :: NTH, TID 94 95! Variables to compute the matrix element of the plane wave 96! (not included in the original vmat subroutine) 97 integer :: irel, iua, irealim, inmp(3) 98 real(dp) :: kxij, aux(2), dist(3), kpoint(3) 99 real(dp) :: dxp(3), displaat(3) 100 real(dp) :: dxpgrid(3,nsp) 101 complex(dp), pointer :: delkmats(:), t_delkmats(:,:) 102 103#ifdef _TRACE_ 104 integer :: MPIerror 105#endif 106 107#ifdef DEBUG 108 call write_debug( ' PRE delk' ) 109#endif 110 111#ifdef _TRACE_ 112 call MPI_Barrier( MPI_Comm_World, MPIerror ) 113 call MPItrace_event( 1000, 4 ) 114#endif 115 116! Start time counter 117 call timer('delk',1) 118 119! Initialize the matrix elements of exp(i*\vec{k} \vec{r}) 120 delkmat(:) = 0.0_dp 121 kpoint(:) = dble(iexpikr) * wavevector(:) 122 123!! For debugging 124! kpoint(:) = 0.0_dp 125!! End debugging 126 127! Find atomic cutoff radii 128 nullify(r2cut) 129 call re_alloc( r2cut, 1, nsmax, 'r2cut', 'delk' ) 130 r2cut = 0.0_dp 131 do i = 1,nuotot 132 ia = iaorb(i) 133 is = isa(ia) 134 io = iphorb(i) 135 r2cut(is) = max( r2cut(is), rcut(is,io)**2 ) 136 end do 137 138! Set algorithm logical 139 ParallelLocal = (Nodes > 1) 140 lenx = meshLim(2,1) - meshLim(1,1) + 1 141 leny = meshLim(2,2) - meshLim(1,2) + 1 142 lenz = meshLim(2,3) - meshLim(1,3) + 1 143 lenxy = lenx*leny 144 145! Find value of maxloc 146 maxloc2 = maxval(endpht(1:np)-endpht(0:np-1)) 147 maxloc = maxloc2 + minloc 148 maxloc = min( maxloc, no ) 149 triang = (maxloc+1)*(maxloc+2)/2 150 if ( ParallelLocal ) then 151 if ( nrowsDscfL > 0 ) then 152 nvmaxl = listdlptr(nrowsDscfL) + numdl(nrowsDscfL) 153 else 154 nvmaxl = 1 155 end if 156 end if 157 158! Allocate local memory 159!$OMP parallel default(shared), & 160!$OMP&shared(NTH,t_DscfL,t_delkmats), & 161!$OMP&private(TID,last,delkmats,irealim), & 162!$OMP&private(ip,nc,nlocal,ic,imp,i,il,iu,iul,ii,ind,j,ijl,jl), & 163!$OMP&private(lasta,lastop,ia,is,iop,isp,ix,dxsp,r2sp,nphiloc,iphi,jc), & 164!$OMP&private(Vij,VClocal,DscfL,ilocal,ilc,iorb,Vlocal,Clocal,phia), & 165!$OMP&private(irel,inmp,dxp,dxpgrid,dist,kxij,iua,displaat,aux) 166 167!$OMP single 168#ifdef _OPENMP 169 NTH = omp_get_num_threads( ) 170#else 171 NTH = 1 172#endif 173!$OMP end single ! implicit barrier, IMPORTANT 174 175#ifdef _OPENMP 176 TID = omp_get_thread_num( ) + 1 177#else 178 TID = 1 179#endif 180 181 nullify(Clocal,phia,ilocal,ilc,iorb,Vlocal) 182!$OMP critical 183! Perhaps the critical section is not needed, 184! however it "tells" the OS to allocate per 185! thread, possibly waiting for each thread to 186! place the memory in the best position. 187 allocate( Clocal(nsp,maxloc2) ) 188 allocate( ilocal(no) , ilc(maxloc2) , iorb(maxloc) ) 189 allocate( Vlocal(triang,2) ) 190 if ( DirectPhi ) allocate( phia(maxoa,nsp) ) 191!$OMP end critical 192 193!$OMP single 194 if ( ParallelLocal ) then 195 nullify( t_DscfL ) 196 call re_alloc( t_DscfL, 1, nvmaxl, 1, 2, 1, NTH, & 197 'DscfL', 'delk' ) 198 else 199 if ( NTH > 1 ) then 200 nullify( t_delkmats ) 201 call re_alloc( t_delkmats, 1, nvmax, 1, NTH, & 202 'delkmats', 'delk' ) 203 end if 204 end if 205!$OMP end single 206 207 if ( ParallelLocal ) then 208 DscfL => t_DscfL(1:nvmaxl,1:2,TID) 209 DscfL(1:nvmaxl,1:2) = 0.0_dp 210 else 211 if ( NTH > 1 ) then 212 delkmats => t_delkmats(1:nvmax,TID) 213 else 214 delkmats => delkmat 215 end if 216 end if 217 218! Full initializations done only once 219 ilocal(1:no) = 0 220 iorb(1:maxloc) = 0 221 Vlocal(1:triang,1:2) = 0.0_dp 222 last = 0 223 224! Loop over grid points 225!$OMP do 226 do ip = 1,np 227! Find number of nonzero orbitals at this point 228 nc = endpht(ip) - endpht(ip-1) 229! Find new required size of Vlocal 230 nlocal = last 231 do ic = 1,nc 232 imp = endpht(ip-1) + ic 233 i = lstpht(imp) 234 if (ilocal(i) .eq. 0) nlocal = nlocal + 1 235 end do 236 237! If overflooded, add Vlocal to delkmat and reinitialize it 238 if (nlocal > maxloc .and. last > 0) then 239 if ( ParallelLocal ) then 240 do il = 1,last 241 i = iorb(il) 242 iu = indxuo(i) 243 iul = NeedDscfL(iu) 244 if ( i == iu ) then 245 do ii = 1, numdl(iul) 246 ind = listdlptr(iul)+ii 247 j = listdl(ind) 248 jl = ilocal(j) 249 ijl = idx_ijl(il,jl) 250! The variables we want to compute in this subroutine are complex numbers 251! Here, when irealim =1 we refer to the real part, and 252! when irealim = 2 we refer to the imaginary part 253 DscfL(ind,:) = DscfL(ind,:) + Vlocal(ijl,:) * dVol 254 end do 255 else 256 ia = iaorb(i) 257 iua = indxua(ia) 258 do ix = 1, 3 259 displaat(ix) = & 260 (iatfold(1,ia)*nmsc(1))*cmesh(ix,1)+ & 261 (iatfold(2,ia)*nmsc(2))*cmesh(ix,2)+ & 262 (iatfold(3,ia)*nmsc(3))*cmesh(ix,3) 263 end do 264 dist(:) = xa(:,iua) - xa(:,ia) - displaat(:) 265 kxij = kpoint(1) * dist(1) + & 266 kpoint(2) * dist(2) + & 267 kpoint(3) * dist(3) 268 do ii = 1, numdl(iul) 269 ind = listdlptr(iul)+ii 270 j = listsc( i, iu, listdl(ind) ) 271 jl = ilocal(j) 272 ijl = idx_ijl(il,jl) 273 aux(1) = Vlocal(ijl,1) * dcos(kxij) - & 274 Vlocal(ijl,2) * dsin(kxij) 275 aux(2) = Vlocal(ijl,2) * dcos(kxij) + & 276 Vlocal(ijl,1) * dsin(kxij) 277 DscfL(ind,:) = DscfL(ind,:) + aux(:) * dVol 278 end do 279 end if 280 end do 281 else 282 283 do il = 1,last 284 i = iorb(il) 285 iu = indxuo(i) 286 call GlobalToLocalOrb( iu, Node, Nodes, iul ) 287 if (i == iu) then 288 do ii = 1, numVs(iul) 289 ind = listVsptr(iul)+ii 290 j = listVs(ind) 291 jl = ilocal(j) 292 ijl = idx_ijl(il,jl) 293 delkmats(ind) = delkmats(ind) + & 294 cmplx(Vlocal(ijl,1), Vlocal(ijl,2), kind=dp) * dVol 295 296 end do 297 else 298 ia = iaorb(i) 299 iua = indxua(ia) 300 do ix = 1, 3 301 displaat(ix) = & 302 (iatfold(1,ia)*nmsc(1))*cmesh(ix,1)+ & 303 (iatfold(2,ia)*nmsc(2))*cmesh(ix,2)+ & 304 (iatfold(3,ia)*nmsc(3))*cmesh(ix,3) 305 end do 306 dist(:) = xa(:,iua) - xa(:,ia) - displaat(:) 307 kxij = kpoint(1) * dist(1) + & 308 kpoint(2) * dist(2) + & 309 kpoint(3) * dist(3) 310 do ii = 1, numVs(iul) 311 ind = listVsptr(iul)+ii 312 j = listsc( i, iu, listVs(ind) ) 313 jl = ilocal(j) 314 ijl = idx_ijl(il,jl) 315 aux(1) = Vlocal(ijl,1) * dcos(kxij) - & 316 Vlocal(ijl,2) * dsin(kxij) 317 aux(2) = Vlocal(ijl,2) * dcos(kxij) + & 318 Vlocal(ijl,1) * dsin(kxij) 319 320 delkmats(ind) = delkmats(ind) + & 321 cmplx(aux(1),aux(2),kind=dp) * dVol 322 323 end do 324 end if 325 end do 326 327 end if 328 329! Reset local arrays 330 do ii = 1, last 331 ilocal(iorb(ii)) = 0 332 end do 333 iorb(1:last) = 0 334 ijl = (last+1)*(last+2)/2 335 Vlocal(1:ijl,1:2) = 0.0_dp 336 last = 0 337 end if 338 339! Look for required orbitals not yet in Vlocal 340 if ( nlocal > last ) then 341 do ic = 1,nc 342 imp = endpht(ip-1) + ic 343 i = lstpht(imp) 344 if (ilocal(i) == 0) then 345 last = last + 1 346 ilocal(i) = last 347 iorb(last) = i 348 end if 349 end do 350 end if 351 352 if ( DirectPhi ) then 353 lasta = 0 354 lastop = 0 355 do ic = 1 , nc 356 imp = endpht(ip-1) + ic 357 i = lstpht(imp) 358 il = ilocal(i) 359 ia = iaorb(i) 360 iop = listp2(imp) 361 ilc(ic) = il 362 363! Localize the position of the mesh point 364 irel = idop(iop) + ipa(ia) 365 call ipack( -1, 3, nem, inmp, irel ) 366 inmp(:) = inmp(:) + (meshLim(1,:) - 1) 367 inmp(:) = inmp(:) - 2 * ne(:) 368 369 dxp(:) = cmesh(:,1) * inmp(1) + & 370 cmesh(:,2) * inmp(2) + & 371 cmesh(:,3) * inmp(3) 372 373 do isp = 1, nsp 374 dxpgrid(:,isp) = dxp(:) + xdsp(:,isp) 375 end do 376 377! Generate phi values 378 if ( ia /= lasta .or. iop /= lastop ) then 379 lasta = ia 380 lastop = iop 381 is = isa(ia) 382 do isp = 1,nsp 383 dxsp(:) = xdsp(:,isp) + xdop(:,iop) - dxa(:,ia) 384 r2sp = sum(dxsp**2) 385 if ( r2sp < r2cut(is) ) then 386!$OMP critical 387 call all_phi( is, +1, dxsp, maxoa, nphiloc, phia(:,isp) ) 388!$OMP end critical 389 else 390 phia(:,isp) = 0.0_dp 391 end if 392 end do 393 end if 394 iphi = iphorb(i) 395 Clocal(:,ic) = phia(iphi,:) 396 397! Pre-multiply V and Clocal(,ic) 398 Vij(:) = 0._dp 399 do isp = 1 , nsp 400 kxij = kpoint(1) * dxpgrid(1,isp) + & 401 kpoint(2) * dxpgrid(2,isp) + & 402 kpoint(3) * dxpgrid(3,isp) 403 VClocal(1,isp) = dcos(kxij) * Clocal(isp,ic) 404 VClocal(2,isp) = dsin(kxij) * Clocal(isp,ic) 405 Vij(:) = Vij(:) + VClocal(:,isp) * Clocal(isp,ic) 406 end do 407 408! ic == jc, hence we simply do 409 ijl = idx_ijl(il,ilc(ic)) 410 Vlocal(ijl,:) = Vlocal(ijl,:) + Vij(:) 411 412! Loop on second orbital of mesh point 413 do jc = 1 , ic - 1 414 jl = ilc(jc) 415 416! Loop over sub-points 417 Vij(:) = 0.0_dp 418 do isp = 1 , nsp 419 Vij(:) = Vij(:) + VClocal(:,isp) * Clocal(isp,jc) 420 end do 421 422 ijl = idx_ijl(il,jl) 423 if ( il == jl ) then 424 Vlocal(ijl,:) = Vlocal(ijl,:) + Vij(:) * 2._dp 425 else 426 Vlocal(ijl,:) = Vlocal(ijl,:) + Vij(:) 427 end if 428 429 end do 430 431 end do 432 433 else 434 435! Loop on first orbital of mesh point 436 do ic = 1 , nc 437 imp = endpht(ip-1) + ic 438 i = lstpht(imp) 439 il = ilocal(i) 440 ia = iaorb(i) 441 iop = listp2(imp) 442 ilc(ic) = il 443 444! Localize the position of the mesh point 445 irel = idop(iop) + ipa(ia) 446 call ipack( -1, 3, nem, inmp, irel ) 447 inmp(:) = inmp(:) + (meshLim(1,:) - 1) 448 inmp(:) = inmp(:) - 2 * ne(:) 449 450 dxp(:) = cmesh(:,1) * inmp(1) + & 451 cmesh(:,2) * inmp(2) + & 452 cmesh(:,3) * inmp(3) 453 454 do isp = 1, nsp 455 dxpgrid(:,isp) = dxp(:) + xdsp(:,isp) 456 end do 457 458! Retrieve phi values 459 Clocal(:,ic) = phi(:,imp) 460 461! Pre-multiply V and Clocal(,ic) 462 Vij(:) = 0._dp 463 do isp = 1 , nsp 464 kxij = kpoint(1) * dxpgrid(1,isp) + & 465 kpoint(2) * dxpgrid(2,isp) + & 466 kpoint(3) * dxpgrid(3,isp) 467 VClocal(1,isp) = dcos(kxij) * Clocal(isp,ic) 468 VClocal(2,isp) = dsin(kxij) * Clocal(isp,ic) 469 Vij(:) = Vij(:) + VClocal(:,isp) * Clocal(isp,ic) 470 end do 471 472! ic == jc, hence we simply do 473 ijl = idx_ijl(il,ilc(ic)) 474 Vlocal(ijl,:) = Vlocal(ijl,:) + Vij(:) 475 476! Loop on second orbital of mesh point 477 do jc = 1 , ic - 1 478 jl = ilc(jc) 479 480! Loop over sub-points 481 Vij(:) = 0.0_dp 482 do isp = 1 , nsp 483 Vij(:) = Vij(:) + VClocal(:,isp) * Clocal(isp,jc) 484 end do 485 486 ijl = idx_ijl(il,jl) 487 if ( il == jl ) then 488 Vlocal(ijl,:) = Vlocal(ijl,:) + Vij(:) * 2._dp 489 else 490 Vlocal(ijl,:) = Vlocal(ijl,:) + Vij(:) 491 end if 492 493 end do 494 495 end do 496 end if 497 498 end do 499!$OMP end do nowait 500 501!$OMP barrier 502 503! Add final Vlocal to delkmat 504 if ( ParallelLocal .and. last > 0 ) then 505 506 do il = 1 , last 507 i = iorb(il) 508 iu = indxuo(i) 509 iul = NeedDscfL(iu) 510 if ( i == iu ) then 511 do ii = 1, numdl(iul) 512 ind = listdlptr(iul)+ii 513 j = listdl(ind) 514 jl = ilocal(j) 515 ijl = idx_ijl(il,jl) 516 DscfL(ind,:) = DscfL(ind,:) + Vlocal(ijl,:) * dVol 517 end do 518 else 519 ia = iaorb(i) 520 iua = indxua(ia) 521 do ix = 1, 3 522 displaat(ix) = & 523 (iatfold(1,ia)*nmsc(1))*cmesh(ix,1)+ & 524 (iatfold(2,ia)*nmsc(2))*cmesh(ix,2)+ & 525 (iatfold(3,ia)*nmsc(3))*cmesh(ix,3) 526 end do 527 dist(:) = xa(:,iua) - xa(:,ia) - displaat(:) 528 kxij = kpoint(1) * dist(1) + & 529 kpoint(2) * dist(2) + & 530 kpoint(3) * dist(3) 531 do ii = 1, numdl(iul) 532 ind = listdlptr(iul)+ii 533 j = listsc( i, iu, listdl(ind) ) 534 jl = ilocal(j) 535 ijl = idx_ijl(il,jl) 536 aux(1) = Vlocal(ijl,1) * dcos(kxij) - & 537 Vlocal(ijl,2) * dsin(kxij) 538 aux(2) = Vlocal(ijl,2) * dcos(kxij) + & 539 Vlocal(ijl,1) * dsin(kxij) 540 DscfL(ind,:) = DscfL(ind,:) + aux(:) * dVol 541 end do 542 end if 543 end do 544 545 else if ( last > 0 ) then 546 547 do il = 1 , last 548 i = iorb(il) 549 iu = indxuo(i) 550 if ( i == iu ) then 551 do ii = 1, numVs(iu) 552 ind = listVsptr(iu)+ii 553 j = listVs(ind) 554 jl = ilocal(j) 555 ijl = idx_ijl(il,jl) 556 delkmats(ind) = delkmats(ind) + & 557 cmplx(Vlocal(ijl,1), Vlocal(ijl,2),kind=dp) * dVol 558 end do 559 else 560 ia = iaorb(i) 561 iua = indxua(ia) 562 do ix = 1, 3 563 displaat(ix) = & 564 (iatfold(1,ia)*nmsc(1))*cmesh(ix,1)+ & 565 (iatfold(2,ia)*nmsc(2))*cmesh(ix,2)+ & 566 (iatfold(3,ia)*nmsc(3))*cmesh(ix,3) 567 end do 568 dist(:) = xa(:,iua) - xa(:,ia) - displaat(:) 569 kxij = kpoint(1) * dist(1) + & 570 kpoint(2) * dist(2) + & 571 kpoint(3) * dist(3) 572 do ii = 1, numVs(iu) 573 ind = listVsptr(iu)+ii 574 j = listsc( i, iu, listVs(ind) ) 575 jl = ilocal(j) 576 ijl = idx_ijl(il,jl) 577 aux(1) = Vlocal(ijl,1) * dcos(kxij) - & 578 Vlocal(ijl,2) * dsin(kxij) 579 aux(2) = Vlocal(ijl,2) * dcos(kxij) + & 580 Vlocal(ijl,1) * dsin(kxij) 581 delkmats(ind) = delkmats(ind) + & 582 cmplx(aux(1),aux(2),kind=dp) * dVol 583 end do 584 end if 585 end do 586 587 end if 588 589!$OMP barrier 590 591 if ( ParallelLocal .and. NTH > 1 ) then 592 do irealim = 1, 2 593!$OMP do 594 do ind = 1, nvmaxl 595 do ii = 2, NTH 596 t_DscfL(ind,irealim,1) = t_DscfL(ind,irealim,1) + & 597 t_DscfL(ind,irealim,ii) 598 end do 599 end do 600!$OMP end do nowait 601 end do 602!$OMP barrier 603 else if ( NTH > 1 ) then 604!$OMP do 605 do ind = 1, nvmax 606 do ii = 1, NTH 607 delkmat(ind) = delkmat(ind) + t_delkmats(ind,ii) 608 end do 609 end do 610!$OMP end do 611 end if 612 613! Free local memory 614 deallocate(Clocal,ilocal,ilc,iorb,Vlocal) 615 if ( DirectPhi ) deallocate( phia ) 616 617!$OMP master 618 if ( ParallelLocal ) then 619! Redistribute Hamiltonian from mesh to orbital based distribution 620 DscfL => t_DscfL(1:nvmaxl,1:2,1) 621 call matrixMtoOC( nvmaxl, nvmax, numVs, listVsptr, nuo, DscfL, delkmat ) 622 call de_alloc( t_DscfL, 'DscfL', 'delk' ) 623 else if ( NTH > 1 ) then 624 call de_alloc( t_delkmats, 'delkmats', 'delk' ) 625 end if 626!$OMP end master 627 628!$OMP end parallel 629 630 call de_alloc( r2cut, 'r2cut', 'delk' ) 631 632#ifdef _TRACE_ 633 call MPI_Barrier( MPI_Comm_World, MPIerror ) 634 call MPItrace_event( 1000, 0 ) 635#endif 636 637 call timer('delk',2) 638 639#ifdef DEBUG 640 call write_debug( ' POS delk' ) 641#endif 642 643 contains 644 645! In any case will the compiler most likely inline this 646! small routine. So it should not pose any problem. 647 pure function idx_ijl(i,j) result(ij) 648 integer, intent(in) :: i,j 649 integer :: ij 650 if ( i > j ) then 651 ij = i * (i + 1)/2 + j + 1 652 else 653 ij = j * (j + 1)/2 + i + 1 654 end if 655 end function idx_ijl 656 657 end subroutine delk 658 659end module m_delk 660