1!=----------------------------------------------------------------------= 2MODULE stick_base 3!=----------------------------------------------------------------------= 4 5 USE fft_param 6 IMPLICIT NONE 7 PRIVATE 8 SAVE 9 10 PUBLIC :: sticks_map_set 11 PUBLIC :: sticks_map_index, sticks_sort_new, sticks_dist_new 12 PUBLIC :: sticks_map, sticks_map_allocate 13 PUBLIC :: sticks_map_deallocate, get_sticks 14 15 TYPE sticks_map 16 LOGICAL :: lgamma=.false. ! if .true. the map has gamma symmetry 17 LOGICAL :: lpara=.false. ! if .true. the map is set for parallel and serial, if .false. only serial 18 INTEGER :: mype=0 ! my task id (starting from 0) 19 INTEGER :: nproc=1 ! number of task (as nproc in fft_type_descriptor) 20 INTEGER :: nyfft=1 ! number processors in y-direction (as nproc2 in fft_type_descriptor) 21 INTEGER, ALLOCATABLE :: iproc(:,:) ! the processor index (as in fft_type_descriptor) 22 INTEGER, ALLOCATABLE :: iproc2(:) ! the Y group processor index (as in fft_type_descriptor) 23#if defined(__MPI) 24 INTEGER :: comm = MPI_COMM_NULL 25#else 26 INTEGER :: comm = 0 ! communicator of the fft gruop 27#endif 28 INTEGER :: nstx=0 ! a safe maximum number of sticks on the map 29 INTEGER :: lb(3)=0 ! map's lower bounds 30 INTEGER :: ub(3)=0 ! map's upper bounds 31 INTEGER, ALLOCATABLE :: idx(:) ! the index of each stick 32 INTEGER, ALLOCATABLE :: ist(:,:) ! the cartesian coordinates of each stick 33 INTEGER, ALLOCATABLE :: stown(:,:) ! the owner of each stick 34 INTEGER, ALLOCATABLE :: indmap(:,:) ! the index of each stick (represented on the map) 35 REAL(DP) :: bg(3,3) ! base vectors, the generators of the mapped space 36 END TYPE 37 38!=----------------------------------------------------------------------= 39CONTAINS 40!=----------------------------------------------------------------------= 41 SUBROUTINE sticks_map_deallocate( smap ) 42 IMPLICIT NONE 43 TYPE( sticks_map ) :: smap 44 !write (6,*) ' inside sticks_map_deallocate'; FLUSH(6) 45 IF( ALLOCATED( smap%iproc ) ) DEALLOCATE( smap%iproc ) 46 IF( ALLOCATED( smap%iproc2 ) ) DEALLOCATE( smap%iproc2 ) 47 IF( ALLOCATED( smap%idx ) ) DEALLOCATE( smap%idx ) 48 IF( ALLOCATED( smap%ist ) ) DEALLOCATE( smap%ist ) 49 IF( ALLOCATED( smap%stown ) ) DEALLOCATE( smap%stown ) 50 IF( ALLOCATED( smap%indmap ) ) DEALLOCATE( smap%indmap ) 51 smap%ub = 0 52 smap%lb = 0 53 smap%nstx = 0 54 END SUBROUTINE sticks_map_deallocate 55 56 SUBROUTINE sticks_map_allocate( smap, lgamma, lpara, nyfft, iproc, iproc2, nr1, nr2, nr3, bg, comm ) 57 ! assign 58 IMPLICIT NONE 59 TYPE( sticks_map ) :: smap 60 LOGICAL, INTENT(IN) :: lgamma 61 LOGICAL, INTENT(IN) :: lpara 62 INTEGER, INTENT(IN) :: nyfft 63 INTEGER, INTENT(IN) :: iproc(:,:) 64 INTEGER, INTENT(IN) :: iproc2(:) 65 INTEGER, INTENT(IN) :: nr1, nr2, nr3 66 INTEGER, INTENT(IN) :: comm 67 REAL(DP), INTENT(IN) :: bg(3,3) 68 INTEGER :: lb(3), ub(3) 69 INTEGER :: nzfft, nstx, ierr 70 INTEGER, ALLOCATABLE :: indmap(:,:), stown(:,:), idx(:), ist(:,:) 71 !write (6,*) ' inside sticks_map_allocate'; FLUSH(6) 72 ub(1) = ( nr1 - 1 ) / 2 73 ub(2) = ( nr2 - 1 ) / 2 74 ub(3) = ( nr3 - 1 ) / 2 75 lb = - ub 76 !write(*,*) 'ub,lb',ub,lb 77 nstx = (ub(1)-lb(1)+1)*(ub(2)-lb(2)+1) ! we stay very large indeed 78 !write(*,*) smap%nstx,nstx 79 IF( smap%nstx == 0 ) THEN 80 ! this map is clean, allocate 81 !write (*,*) ' ! this map is clean, allocate ' 82 smap%mype = 0 83 smap%nproc = 1 84 smap%comm = comm 85#if defined(__MPI) 86 CALL MPI_COMM_RANK( smap%comm, smap%mype, ierr ) 87 CALL MPI_COMM_SIZE( smap%comm, smap%nproc, ierr ) 88#endif 89 smap%lgamma = lgamma 90 smap%lpara = lpara 91 smap%comm = comm 92 smap%nstx = nstx 93 smap%ub = ub 94 smap%lb = lb 95 smap%bg = bg 96 smap%nyfft = nyfft 97 nzfft = smap%nproc / nyfft 98 ALLOCATE ( smap%iproc( nyfft, nzfft ), smap%iproc2 ( smap%nproc ) ) 99 smap%iproc = iproc 100 smap%iproc2= iproc2 101 102 IF( ALLOCATED( smap%indmap ) ) THEN 103 CALL fftx_error__(' sticks_map_allocate ',' indmap already allocated ', 1 ) 104 END IF 105 IF( ALLOCATED( smap%stown ) ) THEN 106 CALL fftx_error__(' sticks_map_allocate ',' stown already allocated ', 1 ) 107 END IF 108 IF( ALLOCATED( smap%idx ) ) THEN 109 CALL fftx_error__(' sticks_map_allocate ',' idx already allocated ', 1 ) 110 END IF 111 IF( ALLOCATED( smap%ist ) ) THEN 112 CALL fftx_error__(' sticks_map_allocate ',' ist already allocated ', 1 ) 113 END IF 114 ALLOCATE( smap%indmap ( lb(1):ub(1), lb(2):ub(2) ) ) 115 ALLOCATE( smap%stown ( lb(1):ub(1), lb(2):ub(2) ) ) 116 ALLOCATE( smap%idx( nstx ) ) 117 ALLOCATE( smap%ist( nstx , 2) ) 118 smap%stown = 0 119 smap%indmap = 0 120 smap%idx = 0 121 smap%ist = 0 122 ELSE IF( smap%nstx < nstx .OR. smap%ub(3) < ub(3) ) THEN 123 ! change the size of the map, but keep the data already there 124 !write (*,*) ' ! change the size of the map, but keep the data already there ' 125 IF( smap%lgamma .neqv. lgamma ) THEN 126 CALL fftx_error__(' sticks_map_allocate ',' changing gamma symmetry not allowed ', 1 ) 127 END IF 128 IF( smap%comm /= comm ) THEN 129 CALL fftx_error__(' sticks_map_allocate ',' changing communicator not allowed ', 1 ) 130 END IF 131 ALLOCATE( indmap ( lb(1):ub(1), lb(2):ub(2) ) ) 132 ALLOCATE( stown ( lb(1):ub(1), lb(2):ub(2) ) ) 133 ALLOCATE( idx( nstx ) ) 134 ALLOCATE( ist( nstx , 2) ) 135 idx = 0 136 ist = 0 137 indmap = 0 138 stown = 0 139 idx( 1:smap%nstx ) = smap%idx 140 ist( 1:smap%nstx, : ) = smap%ist 141 indmap( smap%lb(1):smap%ub(1), smap%lb(2):smap%ub(2) ) = smap%indmap( smap%lb(1):smap%ub(1), smap%lb(2):smap%ub(2) ) 142 stown( smap%lb(1):smap%ub(1), smap%lb(2):smap%ub(2) ) = smap%stown( smap%lb(1):smap%ub(1), smap%lb(2):smap%ub(2) ) 143 DEALLOCATE( smap%indmap ) 144 DEALLOCATE( smap%stown ) 145 DEALLOCATE( smap%idx ) 146 DEALLOCATE( smap%ist ) 147 ALLOCATE( smap%indmap ( lb(1):ub(1), lb(2):ub(2) ) ) 148 ALLOCATE( smap%stown ( lb(1):ub(1), lb(2):ub(2) ) ) 149 ALLOCATE( smap%idx( nstx ) ) 150 ALLOCATE( smap%ist( nstx , 2) ) 151 smap%indmap = indmap 152 smap%stown = stown 153 smap%idx = idx 154 smap%ist = ist 155 DEALLOCATE( indmap ) 156 DEALLOCATE( stown ) 157 DEALLOCATE( idx ) 158 DEALLOCATE( ist ) 159 smap%nstx = nstx 160 smap%ub = ub 161 smap%lb = lb 162 smap%bg = bg 163 smap%nyfft = nyfft 164 smap%iproc = iproc 165 smap%iproc2= iproc2 166 ELSE 167 !write(*,*) 'map is the same ... should we update ub,lb ?' 168 !write(*,*) smap%ub,smap%lb 169 IF( smap%lgamma .neqv. lgamma ) THEN 170 CALL fftx_error__(' sticks_map_allocate ',' changing gamma symmetry not allowed ', 2 ) 171 END IF 172 IF( smap%comm /= comm ) THEN 173 CALL fftx_error__(' sticks_map_allocate ',' changing communicator not allowed ', 1 ) 174 END IF 175 END IF 176 RETURN 177 END SUBROUTINE sticks_map_allocate 178 179 SUBROUTINE sticks_map_set( lgamma, ub, lb, bg, gcut, st, comm ) 180 181 ! .. Compute the basic maps of sticks 182 ! .. st(i,j) will contain the number of G vectors of the stick whose indices are (i,j). 183 184 IMPLICIT NONE 185 186 LOGICAL, INTENT(in) :: lgamma ! if true use gamma point simmetry 187 INTEGER, INTENT(in) :: ub(:) ! upper bounds for i-th grid dimension 188 INTEGER, INTENT(in) :: lb(:) ! lower bounds for i-th grid dimension 189 REAL(DP) , INTENT(in) :: bg(:,:) ! reciprocal space base vectors 190 REAL(DP) , INTENT(in) :: gcut ! cut-off for potentials 191 INTEGER, OPTIONAL, INTENT(in) :: comm ! communicator of the g-vec group 192 ! 193 ! stick map for wave functions, note that map is taken in YZ plane 194 ! 195 INTEGER, INTENT(out) :: st( lb(1) : ub(1), lb(2) : ub(2) ) 196 REAL(DP) :: b1(3), b2(3), b3(3) 197 INTEGER :: i1, i2, i3, n1, n2, n3, mype, nproc, ierr 198 REAL(DP) :: amod 199 INTEGER :: ngm 200 !write (6,*) 'inside sticks_map_set gcut=',gcut; FLUSH(6) 201 !write (6,*) ub,lb 202 203 st = 0 204 b1(:) = bg(:,1) 205 b2(:) = bg(:,2) 206 b3(:) = bg(:,3) 207 208 n1 = max( abs( lb(1) ), abs( ub(1) ) ) 209 n2 = max( abs( lb(2) ), abs( ub(2) ) ) 210 n3 = max( abs( lb(3) ), abs( ub(3) ) ) 211 212 mype = 0 213 nproc = 1 214#if defined(__MPI) 215 IF( PRESENT( comm ) ) THEN 216 CALL MPI_COMM_RANK( comm, mype, ierr ) 217 CALL MPI_COMM_SIZE( comm, nproc, ierr ) 218 END IF 219#endif 220 !write (6,*) ' N1,N2,N3 ', n1,n2,n3 221 ngm=0 222 loop1: DO i1 = - n1, n1 223 ! 224 ! Gamma-only: exclude space with x<0 225 ! 226 IF ( (lgamma .and. i1 < 0) .OR. ( MOD( i1 + n1, nproc ) /= mype )) CYCLE loop1 227 ! 228 loop2: DO i2 = - n2, n2 229 ! 230 ! Gamma-only: exclude plane with x=0, y<0 231 ! 232 IF(lgamma .and. i1 == 0.and. i2 < 0) CYCLE loop2 233 ! 234 loop3: DO i3 = - n3, n3 235 ! 236 ! Gamma-only: exclude line with x=0, y=0, z<0 237 ! 238 IF(lgamma .and. i1 == 0 .and. i2 == 0 .and. i3 < 0) CYCLE loop3 239 ! 240 amod = (i1 * b1 (1) + i2 * b2 (1) + i3 * b3 (1) ) **2 + & 241 (i1 * b1 (2) + i2 * b2 (2) + i3 * b3 (2) ) **2 + & 242 (i1 * b1 (3) + i2 * b2 (3) + i3 * b3 (3) ) **2 243 IF (amod <= gcut ) then 244 st( i1, i2 ) = st( i1, i2 ) + 1 245 ngm =ngm + 1 246 endif 247 ENDDO loop3 248 ENDDO loop2 249 ENDDO loop1 250 251#if defined(__MPI) 252 IF( PRESENT( comm ) ) THEN 253 CALL MPI_ALLREDUCE(MPI_IN_PLACE, st, size(st), MPI_INTEGER, MPI_SUM, comm, ierr) 254 CALL MPI_ALLREDUCE(MPI_IN_PLACE, ngm, 1, MPI_INTEGER, MPI_SUM, comm, ierr) 255 END IF 256#endif 257 258 !write (6,*) ' NGM ', ngm 259 260 261 RETURN 262 END SUBROUTINE sticks_map_set 263 264!=----------------------------------------------------------------------= 265 266 SUBROUTINE sticks_map_index( ub, lb, st, in1, in2, ngc, index_map ) 267 268 IMPLICIT NONE 269 270 INTEGER, INTENT(in) :: ub(:), lb(:) 271 INTEGER, INTENT(in) :: st( lb(1): ub(1), lb(2):ub(2) ) ! stick map for potential 272 INTEGER, INTENT(inout) :: index_map( lb(1): ub(1), lb(2):ub(2) ) ! keep track of sticks index 273 274 INTEGER, INTENT(out) :: in1(:), in2(:) 275 INTEGER, INTENT(out) :: ngc(:) 276 277 INTEGER :: j1, j2, i1, i2, i3, nct, min_size, ind 278 LOGICAL :: ok 279 !write (6,*) ' inside sticks_map_index'; FLUSH(6) 280 281! 282! ... initialize the sticks indexes array ist 283! ... nct counts columns containing G-vectors for the dense grid 284! ... ncts counts columns contaning G-vectors for the smooth grid 285! 286 nct = MAXVAL( index_map ) 287 ngc = 0 288 289 min_size = min( size( in1 ), size( in2 ), size( ngc ) ) 290 291 DO j2 = 0, ( ub(2) - lb(2) ) 292 DO j1 = 0, ( ub(1) - lb(1) ) 293 i1 = j1 294 IF( i1 > ub(1) ) i1 = lb(1) + ( i1 - ub(1) ) - 1 295 i2 = j2 296 IF( i2 > ub(2) ) i2 = lb(2) + ( i2 - ub(2) ) - 1 297 IF( st( i1, i2 ) > 0 ) THEN 298 IF( index_map( i1, i2 ) == 0 ) THEN 299 nct = nct + 1 300 index_map( i1, i2 ) = nct 301 END IF 302 ind = index_map( i1, i2 ) 303 IF( nct > min_size ) & 304 CALL fftx_error__(' sticks_map_index ',' too many sticks ', nct ) 305 in1(ind) = i1 306 in2(ind) = i2 307 ngc(ind) = st( i1 , i2) 308 ENDIF 309 ENDDO 310 ENDDO 311 312 RETURN 313 END SUBROUTINE sticks_map_index 314 315!=----------------------------------------------------------------------= 316 317 SUBROUTINE sticks_sort_new( parallel, ng, nct, idx ) 318 319! ... This subroutine sorts the sticks indexes, according to 320! ... the length and type of the sticks, wave functions sticks 321! ... first, then smooth mesh sticks, and finally potential 322! ... sticks 323 324 ! lengths of sticks, ngc for potential mesh, ngcw for wave functions mesh 325 ! and ngcs for smooth mesh 326 327 IMPLICIT NONE 328 329 LOGICAL, INTENT(in) :: parallel 330 INTEGER, INTENT(in) :: ng(:) 331 332 ! nct, total number of sticks 333 334 INTEGER, INTENT(in) :: nct 335 336 ! index, on output, new sticks indexes 337 338 INTEGER, INTENT(inout) :: idx(:) 339 340 INTEGER :: mc, ic, nc 341 INTEGER, ALLOCATABLE :: iaux(:) 342 INTEGER, ALLOCATABLE :: itmp(:) 343 REAL(DP), ALLOCATABLE :: aux(:) 344 !write (6,*) ' inside sticks_map_sort_new'; FLUSH(6) 345 346 ! we need to avoid sorting elements already sorted previously 347 ! build inverse indexes 348 ALLOCATE( iaux( nct ) ) 349 iaux = 0 350 DO mc = 1, nct 351 IF( idx( mc ) > 0 ) iaux( idx( mc ) ) = mc 352 END DO 353 ! 354 ! check idx has no "hole" 355 ! 356 IF( idx( 1 ) == 0 ) THEN 357 ic = 0 358 DO mc = 2, nct 359 IF( idx( mc ) /= 0 ) THEN 360 CALL fftx_error__(' sticks_sort ',' non contiguous indexes 1 ', nct ) 361 END IF 362 END DO 363 ELSE 364 ic = 1 365 DO mc = 2, nct 366 IF( idx( mc ) == 0 ) EXIT 367 ic = ic + 1 368 END DO 369 DO mc = ic+1, nct 370 IF( idx( mc ) /= 0 ) THEN 371 CALL fftx_error__(' sticks_sort ',' non contiguous indexes 2 ', nct ) 372 END IF 373 END DO 374 END IF 375 376 IF( parallel ) THEN 377 ALLOCATE( aux( nct ) ) 378 ALLOCATE( itmp( nct ) ) 379 itmp = 0 380 nc = 0 381 DO mc = 1, nct 382 IF( ng( mc ) > 0 .AND. iaux( mc ) == 0 ) THEN 383 nc = nc + 1 384 aux( nc ) = -ng(mc) 385 itmp( nc ) = mc 386 END IF 387 ENDDO 388 CALL hpsort( nc, aux, itmp) 389 DO mc = 1, nc 390 idx( ic + mc ) = itmp( mc ) 391 END DO 392 DEALLOCATE( itmp ) 393 DEALLOCATE( aux ) 394 ELSE 395 DO mc = 1, nct 396 IF( ng(mc) > 0 .AND. iaux(mc) == 0 ) THEN 397 ic = ic + 1 398 idx(ic) = mc 399 ENDIF 400 ENDDO 401 ENDIF 402 403 DEALLOCATE( iaux ) 404 405 RETURN 406 END SUBROUTINE sticks_sort_new 407 408!=----------------------------------------------------------------------= 409 410 SUBROUTINE sticks_dist_new( lgamma, mype, nproc, nyfft, iproc, iproc2, ub, lb, idx, in1, in2, ngc, nct, ncp, ngp, stown, ng ) 411 412 IMPLICIT NONE 413 414 LOGICAL, INTENT(in) :: lgamma 415 INTEGER, INTENT(in) :: mype 416 INTEGER, INTENT(in) :: nproc 417 INTEGER, INTENT(in) :: nyfft 418 INTEGER, INTENT(in) :: iproc(:,:),iproc2(:) 419 420 INTEGER, INTENT(in) :: ub(:), lb(:), idx(:) ! map's limits 421 INTEGER, INTENT(inout) :: stown(lb(1): ub(1), lb(2):ub(2) ) ! stick owner's map, including previously assigned ones 422 INTEGER, INTENT(in) :: in1(:), in2(:) ! cartesian coordinate of each column, ordered according to idx index 423 INTEGER, INTENT(in) :: ngc(:) ! number of G-vector of each column, ordered according to idx index 424 INTEGER, INTENT(in) :: nct ! total number of relevant sticks 425 INTEGER, INTENT(out):: ncp(:) ! number of sticks (columns) per processor 426 INTEGER, INTENT(out):: ngp(:) ! number of G-vectors per processor 427 INTEGER, INTENT(out):: ng ! number of G-vector of this processor. that is ng = ngp(mype+1) 428 429 INTEGER :: mc, i1, i2, j, jj, icnt, gr, j2, j3 430 INTEGER, ALLOCATABLE :: yc(:), yg(:) ! number of relevant columns and G-vectors in a given yz plane 431 INTEGER, ALLOCATABLE :: ygr(:) ! element in the nyfft group to which a YZ plane belong 432 INTEGER, ALLOCATABLE :: ygrp(:), ygrc(:), ygrg(:) ! number of yz planes, relevant columns and G-vectors belonging to a given element in the nyfft group 433 !write (6,*) ' inside sticks_map_dist_new'; FLUSH(6) 434 435 ! distribute X values first 436 !write (6,*) 'nyfft inside sticks_dist_new', nyfft 437 ALLOCATE ( yc(lb(1):ub(1)), yg(lb(1):ub(1)), ygr(lb(1):ub(1)) ) ; yc = 0; yg = 0 ; ygr = 0 438 ALLOCATE ( ygrp(nyfft), ygrc(nyfft), ygrg(nyfft) ) ; ygrp = 0; ygrc = 0; ygrg = 0 439 DO mc =1, nct 440 if( idx( mc ) < 1 ) CYCLE 441 i1 = in1( idx( mc ) ) 442 i2 = in2( idx( mc ) ) 443 IF ( ngc( idx(mc) ) > 0 ) THEN 444 yc(i1) = yc(i1) + 1 445 yg(i1) = yg(i1) + ngc(idx(mc)) 446 END IF 447 IF (stown(i1,i2) > 0) THEN 448 gr = iproc2(stown(i1,i2)) 449 IF (ygr(i1) == 0 ) ygr(i1) = gr 450 if (ygr(i1) .ne. gr ) CALL fftx_error__(' sticks_dist ',' ygroups are not compatible ', 1 ) 451 END IF 452 END DO 453 do i1 = lb(1),ub(1) 454 if (ygr(i1) == 0 ) cycle 455 ygrp(ygr(i1)) = ygrp(ygr(i1)) + 1 456 ygrc(ygr(i1)) = ygrc(ygr(i1)) + yc(i1) 457 ygrg(ygr(i1)) = ygrg(ygr(i1)) + yg(i1) 458 end do 459 460 !write (6,*) 'inside sticks_dist_new', lb(1),ub(1) 461 !write (6,'(4i4)') (i1,yc(i1),yg(i1),ygr(i1),i1=lb(1),ub(1)) 462 !write (6,*) 'initial nyfft distribution' 463 !write (6,'(4i5)') (i1,ygrp(i1),ygrc(i1),ygrg(i1),i1=1,nyfft) 464 465 ncp = 0 466 ngp = 0 467 icnt = 0 468 469 DO mc = 1, nct 470 471 IF( idx( mc ) < 1 ) CYCLE 472 473 i1 = in1( idx( mc ) ) 474 i2 = in2( idx( mc ) ) 475! 476 IF ( lgamma .and. ( (i1 < 0) .or. ( (i1 == 0) .and. (i2 < 0) ) ) ) GOTO 30 477! 478 IF (ygr(i1) == 0 ) THEN 479 j2=1 480 DO j=1,nyfft 481 IF ( ygrg(j) < ygrg(j2) ) THEN 482 j2 = j 483 ELSEIF ( ( ygrg(j) == ygrg(j2) ) .and. ( ygrc(j) < ygrc(j2) ) ) THEN 484 j2 = j 485 END IF 486 END DO 487 ygr(i1) = j2 488 ygrp(j2) = ygrp(j2) + 1 489 ygrc(j2) = ygrc(j2) + yc(i1) 490 ygrg(j2) = ygrg(j2) + yg(i1) 491 ELSE 492 j2=ygr(i1) 493 END IF 494 495 IF ( ngc( idx(mc) ) > 0 .AND. stown(i1,i2) == 0 ) THEN 496 jj = iproc(j2,1) 497 DO j3 = 1, nproc / nyfft 498 j = iproc(j2,j3) 499 IF ( ngp(j) < ngp(jj) ) THEN 500 jj = j 501 ELSEIF ( ( ngp(j) == ngp(jj) ) .and. ( ncp(j) < ncp(jj) ) ) THEN 502 jj = j 503 ENDIF 504 ENDDO 505 stown(i1,i2) = jj 506 END IF 507 IF ( ngc( idx(mc) ) > 0 ) THEN 508 ncp( stown(i1,i2) ) = ncp( stown(i1,i2) ) + 1 509 ngp( stown(i1,i2) ) = ngp( stown(i1,i2) ) + ngc( idx(mc) ) 510 ENDIF 511 30 CONTINUE 512 ENDDO 513 ! 514 ng = ngp( mype + 1 ) 515 ! 516 IF ( lgamma ) THEN 517 ! when gamma symmetry is used only the sticks of half reciprocal space 518 ! are generated, then here we pair-up the sticks with those of the other 519 ! half of the space, using the gamma symmetry relation 520 ! Note that the total numero of stick "nct" is not modified 521 DO mc = 1, nct 522 IF( idx( mc ) < 1 ) CYCLE 523 IF( ngc( idx(mc) ) < 1 ) CYCLE 524 i1 = in1( idx(mc) ) 525 i2 = in2( idx(mc) ) 526 IF( i1 == 0 .and. i2 == 0 ) THEN 527 jj = stown( i1, i2 ) 528 IF( jj > 0 ) ngp( jj ) = ngp( jj ) + ngc( idx(mc) ) - 1 529 ELSE 530 jj = stown( i1, i2 ) 531 IF( jj > 0 ) THEN 532 stown( -i1, -i2 ) = jj 533 ncp( jj ) = ncp( jj ) + 1 534 ngp( jj ) = ngp( jj ) + ngc( idx(mc) ) 535 ENDIF 536 ENDIF 537 ENDDO 538 ENDIF 539 540 !write (6,*) 'final nyfft distribution' 541 !write (6,'(4i5)') (i1,ygrp(i1),ygrc(i1),ygrg(i1),i1=1,nyfft) 542 DEALLOCATE ( ygrp, ygrc, ygrg ) 543 DEALLOCATE ( yc, yg, ygr ) 544 545 RETURN 546 END SUBROUTINE sticks_dist_new 547 548!--------------------------------------------------------------------- 549 550 SUBROUTINE get_sticks( smap, gcut, nstp, sstp, st, nst, ng ) 551 552 IMPLICIT NONE 553 554 TYPE( sticks_map ), INTENT(INOUT) :: smap 555 REAL(DP) , INTENT(in) :: gcut ! kinetic energy cut-off for this stick distribution 556 557 INTEGER, INTENT(out) :: st(smap%lb(1): smap%ub(1), smap%lb(2):smap%ub(2) ) 558 ! in output it contains (if>0) the processor_id+1 of the owner of a given stick 559 ! internally it's used to contain the number of G-vectors of a given stick 560 INTEGER, INTENT(out) :: nstp(:) ! number of sticks per processor 561 INTEGER, INTENT(out) :: sstp(:) ! number of G-vectors per processor 562 INTEGER, INTENT(out) :: nst ! total number of sticks 563 INTEGER, INTENT(out) :: ng ! number of G-vector of this processor. that is ng = sstp(mype+1) 564 565 INTEGER, ALLOCATABLE :: ngc(:) 566 INTEGER :: ic 567 !write (6,*) ' inside get_sticks gcut=',gcut; FLUSH(6) 568 !write (6,*) smap%lb, smap%ub 569 570 IF( .NOT. ALLOCATED( smap%stown ) ) THEN 571 CALL fftx_error__(' get_sticks ',' sticks map, not allocated ', 1 ) 572 END IF 573 574 st = 0 575 CALL sticks_map_set( smap%lgamma, smap%ub, smap%lb, smap%bg, gcut, st, smap%comm ) 576 577 ALLOCATE( ngc ( SIZE( smap%idx ) ) ) 578 ngc = 0 579 CALL sticks_map_index( smap%ub, smap%lb, st, smap%ist(:,1), smap%ist(:,2), ngc, smap%indmap ) 580 nst = count( st > 0 ) 581 CALL sticks_sort_new( smap%nproc>1, ngc, SIZE(smap%idx), smap%idx ) 582 583 CALL sticks_dist_new( smap%lgamma, smap%mype, smap%nproc, smap%nyfft, smap%iproc, smap%iproc2, & 584 smap%ub, smap%lb, smap%idx, & 585 smap%ist(:,1), smap%ist(:,2), ngc, SIZE(smap%idx), nstp, sstp, smap%stown, ng ) 586 587 ! assign the owner of each (relavant) stick 588 st = 0 589 DO ic = 1, SIZE( smap%idx ) 590 IF( smap%idx( ic ) > 0 ) THEN 591 IF( ngc( smap%idx( ic ) ) > 0 ) THEN 592 st( smap%ist(smap%idx( ic ),1), smap%ist(smap%idx( ic ),2) ) = & 593 smap%stown( smap%ist(smap%idx( ic ),1),smap%ist(smap%idx( ic ),2)) 594 IF(smap%lgamma) st(-smap%ist(smap%idx( ic),1),-smap%ist(smap%idx( ic ),2)) = & 595 smap%stown( smap%ist(smap%idx( ic ),1),smap%ist(smap%idx( ic ),2)) 596 END IF 597 END IF 598 END DO 599 DEALLOCATE( ngc ) 600 RETURN 601 END SUBROUTINE get_sticks 602 603!--------------------------------------------------------------------- 604 SUBROUTINE hpsort (n, ra, ind) 605 !--------------------------------------------------------------------- 606 ! sort an array ra(1:n) into ascending order using heapsort algorithm. 607 ! n is input, ra is replaced on output by its sorted rearrangement. 608 ! create an index table (ind) by making an exchange in the index array 609 ! whenever an exchange is made on the sorted data array (ra). 610 ! in case of equal values in the data array (ra) the values in the 611 ! index array (ind) are used to order the entries. 612 ! if on input ind(1) = 0 then indices are initialized in the routine, 613 ! if on input ind(1) != 0 then indices are assumed to have been 614 ! initialized before entering the routine and these 615 ! indices are carried around during the sorting process 616 ! 617 ! no work space needed ! 618 ! free us from machine-dependent sorting-routines ! 619 ! 620 ! adapted from Numerical Recipes pg. 329 (new edition) 621 ! 622 IMPLICIT NONE 623 !-input/output variables 624 INTEGER :: n 625 INTEGER :: ind (n) 626 REAL(DP) :: ra (n) 627 !-local variables 628 INTEGER :: i, ir, j, l, iind 629 REAL(DP) :: rra 630 ! 631 IF (n < 1 ) RETURN 632 ! initialize index array 633 IF (ind (1) ==0) THEN 634 DO i = 1, n 635 ind (i) = i 636 ENDDO 637 ENDIF 638 ! nothing to order 639 IF (n < 2) RETURN 640 ! initialize indices for hiring and retirement-promotion phase 641 l = n / 2 + 1 642 ir = n 64310 CONTINUE 644 ! still in hiring phase 645 IF (l>1) THEN 646 l = l - 1 647 rra = ra (l) 648 iind = ind (l) 649 ! in retirement-promotion phase. 650 ELSE 651 ! clear a space at the end of the array 652 rra = ra (ir) 653 ! 654 iind = ind (ir) 655 ! retire the top of the heap into it 656 ra (ir) = ra (1) 657 ! 658 ind (ir) = ind (1) 659 ! decrease the size of the corporation 660 ir = ir - 1 661 ! done with the last promotion 662 IF (ir==1) THEN 663 ! the least competent worker at all ! 664 ra (1) = rra 665 ! 666 ind (1) = iind 667 RETURN 668 ENDIF 669 ENDIF 670 ! wheter in hiring or promotion phase, we 671 i = l 672 ! set up to place rra in its proper level 673 j = l + l 674 ! 675 DO WHILE (j<=ir) 676 IF (j<ir) THEN 677 ! compare to better underling 678 IF (ra (j) <ra (j + 1) ) THEN 679 j = j + 1 680 ELSEIF (ra (j) ==ra (j + 1) ) THEN 681 IF (ind (j) <ind (j + 1) ) j = j + 1 682 ENDIF 683 ENDIF 684 ! demote rra 685 IF (rra<ra (j) ) THEN 686 ra (i) = ra (j) 687 ind (i) = ind (j) 688 i = j 689 j = j + j 690 ELSEIF (rra==ra (j) ) THEN 691 ! demote rra 692 IF (iind<ind (j) ) THEN 693 ra (i) = ra (j) 694 ind (i) = ind (j) 695 i = j 696 j = j + j 697 ELSE 698 ! set j to terminate do-while loop 699 j = ir + 1 700 ENDIF 701 ! this is the right place for rra 702 ELSE 703 ! set j to terminate do-while loop 704 j = ir + 1 705 ENDIF 706 ENDDO 707 ra (i) = rra 708 ind (i) = iind 709 GOTO 10 710 ! 711 END SUBROUTINE hpsort 712 713!=----------------------------------------------------------------------= 714END MODULE stick_base 715!=----------------------------------------------------------------------= 716