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! 8! This code segment has been fully created by: 9! Nick Papior Andersen, 2013, nickpapior@gmail.com 10! Please conctact the author, prior to re-using this code. 11 12! These routines should ALWAYS be kept separable from the SIESTA package. 13! No routines must ever be dependent on any SIESTA MODULES. 14 15! A module that supplements routines that generally are missing in 16! The FORTRAN standard. 17! Some are also more "exotic" in terms of their use and the FORTRAN standard. 18 19! DEVELOPER REMARK: 20! Many of these functions are not GNU compliant for nested calls. 21! I.e. SORT(UNIQ(arr)) is NOT allowed and will result in errorneous 22! results. 23 24! The following routines are supplied: 25! - VNORM: perform euclidean norm calculations on vectors or matrices. 26! for complex numbers it is equilvalent to VNORM(abs(z)) 27! - SORT : allows sorting of an array. It currently only exists 28! for an integer array. However, this is also the most stringent 29! case as there has not to be any margin of error (EPS) 30! - SORT_QUICK : allows sorting of an array (using the quick-sort algorithm) 31! - INDEX_SORT_HEAP : allows returning an index array that sorts the input array (using the heap-sort algorithm) 32! - UNIQ : returns the unique integer values of an array 33! - UNIQC: returns the number of unique integer vaules of an array 34! - SFIND: will return the index of an integer in a sorted array (by logical 35! fast stepping). It has optional parameters to search: 36! * NEAREST=+1: if not found, returns the first index HIGHER 37! than the sought value (first index if first 38! value is larger than the searched value) 39! * NEAREST=-1: if not found, returns the first index LOWER 40! than the sought value (last index if last 41! value is larger than the searched value) 42! * NEAREST= 0: if not found, returns -1 (NORMAL behaviour) 43! - MODP : Returns "FORTRAN" indexed MOD, i.e. if MOD would return 0, 44! MODP returns the divisor. 45! - MODULOP : Returns "FORTRAN" indexed MODULO, i.e. if MODULO would return 0, 46! MODULOP returns the divisor. 47! - EYE : a routine to generate identity matrices 48! this is a subroutine as I think a function would produce 49! a large memory foot-print before it actually saves to the array. 50! When called on a complex array, only the REAL diagonal is 1. 51! - TRANSPOSE : transpose matrices 52! - TRACE : transpose matrices 53! - ROTATE : returns a vector rotated theta angles (in radians). 54! For 2D points this rotation is implicitly directioned (in the plane) 55! For 3D points a direction is [1,2,3] is needed to specify the rotation 56! direction. 57! - PROJ : A projection method to project a vector onto a certain subspace 58! Exists only for 3D spaces. 59! - VEC_PROJ : A projection method to project a vector onto another vector 60 61module intrinsic_missing 62 63 implicit none 64 65 private 66 67 integer, parameter :: sp = selected_real_kind(5,10) 68 integer, parameter :: dp = selected_real_kind(10,100) 69 70 ! Naming scheme must be "different" than what could be supplied by 71 ! the standard, for instance, NORM2 is part of F08STD! 72 ! VNORM calculates the euclidiean norm of a vector or matrix 73 interface VNORM 74 module procedure VNORM_sp_1 ,VNORM_sp_2 75 module procedure VNORM_dp_1 ,VNORM_dp_2 76 module procedure VNORM_cp_1 ,VNORM_cp_2 77 module procedure VNORM_zp_1 ,VNORM_zp_2 78 end interface 79 80! Pure functions.. (i.e. callable in interface declarations..) 81 public :: VNORM 82 public :: FSORT, SORT, SORT_QUICK 83 public :: INDEX_SORT_HEAP 84 interface INDEX_SORT_HEAP 85 module procedure INDEX_SORT_HEAP_I4 86 end interface INDEX_SORT_HEAP 87 public :: UNIQ, UNIQC 88 public :: SFIND 89 interface SFIND 90 module procedure SORTED_FIND, SORTED_CLOSE 91 end interface SFIND 92#ifdef INTRINSIC_MISSING_TEST 93 public :: SORTED_FIND 94 public :: SORTED_REC_FIND 95 public :: SORTED_BINARY_FIND 96#endif 97 98 99! Elemental functions (can be called on arrays) 100 public :: MODP, MODULOP 101 102! Missing matrix stuff 103 public :: EYE 104 interface EYE 105 module procedure EYE_i_2D 106 module procedure EYE_sp_2D, EYE_dp_2D 107 module procedure EYE_cp_2D, EYE_zp_2D 108 module procedure EYE_i_1D 109 module procedure EYE_sp_1D, EYE_dp_1D 110 module procedure EYE_cp_1D, EYE_zp_1D 111 end interface 112 113 public :: TRANSPOSE 114 interface TRANSPOSE 115 module procedure TRANSPOSE_dp_2D 116 module procedure TRANSPOSE_zp_2D 117 module procedure TRANSPOSE_dp_1D 118 module procedure TRANSPOSE_zp_1D 119 end interface 120 121 public :: TRACE 122 interface TRACE 123 module procedure TRACE_dp_2D 124 module procedure TRACE_zp_2D 125 module procedure TRACE_dp_1D 126 module procedure TRACE_zp_1D 127 end interface 128 129 ! Calculate the angle between two vectors 130 public :: ANGLE 131 interface ANGLE 132 module procedure ANGLE_sp 133 module procedure ANGLE_dp 134 end interface 135 136! Rotation of points 137 public :: ROTATE 138 interface ROTATE 139 module procedure ROTATE_2D 140 module procedure ROTATE_3D 141 end interface 142 143 ! Projection of N-D vector to N-D space 144 public :: SPC_PROJ 145 interface SPC_PROJ 146 module procedure SPC_PROJ_sp 147 module procedure SPC_PROJ_dp 148 end interface 149 150 ! Index of projection of N-D vector to N-D space 151 public :: IDX_SPC_PROJ 152 interface IDX_SPC_PROJ 153 module procedure IDX_SPC_PROJ_sp 154 module procedure IDX_SPC_PROJ_dp 155 end interface 156 157 ! Projection of a N-D vector onto N-D vector 158 public :: VEC_PROJ 159 interface VEC_PROJ 160 module procedure VEC_PROJ_sp 161 module procedure VEC_PROJ_dp 162 end interface 163 164 ! Projection of a N-D vector onto N-D vector 165 ! by keeping the scaling of the onto projected vector 166 public :: VEC_PROJ_SCA 167 interface VEC_PROJ_SCA 168 module procedure VEC_PROJ_SCA_sp 169 module procedure VEC_PROJ_SCA_dp 170 end interface 171 172contains 173 174 175! ROTATE a point around origo with some angle \theta 176! This is a purely plane rotation 177 pure subroutine ROTATE_2D(v,theta) 178 real(dp), intent(inout) :: v(2) 179 real(dp), intent(in) :: theta 180 real(dp) :: rmT(2,2), vv(2) 181 rmT(1,1) = cos(theta) 182 rmT(1,2) = sin(theta) 183 rmT(2,1) = -rmT(1,2) 184 rmT(2,2) = rmT(1,1) 185 186 vv(1) = sum(rmT(:,1) * v) 187 vv(2) = sum(rmT(:,2) * v) 188 v = vv 189 190 end subroutine ROTATE_2D 191 192 pure subroutine ROTATE_3D(v,theta,dir) 193 real(dp), intent(inout) :: v(3) 194 real(dp), intent(in) :: theta 195 integer, intent(in) :: dir 196 real(dp) :: rmT(3,3), vv(3) 197 rmT(:,:) = 0._dp 198 if ( dir == 3 ) then 199 rmT(1,1) = cos(theta) 200 rmT(1,2) = sin(theta) 201 rmT(2,1) = -rmT(1,2) 202 rmT(2,2) = rmT(1,1) 203 rmT(3,3) = 1._dp 204 else if ( dir == 2 ) then 205 rmT(1,1) = cos(theta) 206 rmT(1,3) = sin(theta) 207 rmT(2,2) = 1._dp 208 rmT(3,1) = -rmT(1,3) 209 rmT(3,3) = rmT(1,1) 210 else if ( dir == 1 ) then 211 rmT(1,1) = 1._dp 212 rmT(2,2) = cos(theta) 213 rmT(2,3) = sin(theta) 214 rmT(3,2) = -rmT(2,3) 215 rmT(3,3) = rmT(2,2) 216 end if 217 218 vv(1) = sum(rmT(:,1) * v) 219 vv(2) = sum(rmT(:,2) * v) 220 vv(3) = sum(rmT(:,3) * v) 221 v = vv 222 223 end subroutine ROTATE_3D 224 225 226! A MOD function which behaves differently on the edge. 227! It will NEVER return 0, but instead return the divisor. 228! This makes it useful in FORTRAN do-loops which are not zero-based 229! but 1-based. 230! Thus we have the following scheme: 231! MODP(x',x) == x for x' == x 232! MODP(x',x) == MOD(x',x) for x' /= x 233 elemental function MODP(a,p) 234 integer, intent(in) :: a,p 235 integer :: MODP 236 MODP = MOD(a-1,p) + 1 237 end function MODP 238 239! A MODULO function which behaves differently on the edge. 240! It will NEVER return 0, but instead return the divisor. 241! This makes it useful in FORTRAN do-loops which are not zero-based 242! but 1-based. 243! Thus we have the following scheme: 244! MODULOP(x',x) == x for x' == x 245! MODULOP(x',x) == MODULOP(x',x) for x' /= x 246 elemental function MODULOP(a,p) 247 integer, intent(in) :: a,p 248 integer :: MODULOP 249 MODULOP = MODULO(a-1,p) + 1 250 end function MODULOP 251 252! Function to return the unique COUNT of an integer array. 253! Thus will return how many DIFFERENT entries there exists. 254 pure function UNIQC(array) 255 ! We dont need to see the size of the array 256 integer, intent(in) :: array(:) 257 integer :: dummy(ubound(array,dim=1)) 258 integer :: UNIQC 259 integer :: i, DA 260 261 DA = ubound(array,dim=1) 262 263 ! The first entry is of course not found anywhere else 264 if ( DA == 0 ) then 265 UNIQC = 0 266 return 267 end if 268 269 ! Everything else has to be checked... 270 UNIQC = 1 271 ! Save the first entry 272 dummy(1) = array(1) 273 do i = 2 , DA 274 ! Apparently, inserting a do-loop is worse than count? 275 ! This makes me worried.... ??? 276! do ii = 1 , UNIQC 277! if ( dummy(ii) == array(i) ) cycle find_num 278! end do 279! UNIQC = UNIQC + 1 280! dummy(UNIQC) = array(i) 281 282 if ( count(dummy(1:UNIQC) == array(i)) == 0 ) then 283 UNIQC = UNIQC + 1 284 dummy(UNIQC) = array(i) 285 end if 286 end do 287 288 end function UNIQC 289 290! Function to return the unique COUNT of an integer array (THIS ARRAY HAS 291! TO BE SORTED). 292! Thus will return how many DIFFERENT entries there exists. 293 pure function UNIQC_SORTED(array) result(UNIQC) 294 ! We dont need to see the size of the array 295 integer, intent(in) :: array(:) 296 integer :: UNIQC 297 integer :: i, DA 298 299 DA = ubound(array,dim=1) 300 301 ! The first entry is of course not found anywhere else 302 if ( DA == 0 ) then 303 UNIQC = 0 304 return 305 end if 306 307 ! Everything else has to be checked... 308 UNIQC = 1 309 ! Save the first entry 310 do i = 2 , DA 311 if ( array(i-1) /= array(i) ) then 312 UNIQC = UNIQC + 1 313 end if 314 end do 315 316 end function UNIQC_SORTED 317 318 pure function UNIQ(array) 319 integer, intent(in) :: array(:) 320 integer :: UNIQ(UNIQC(array)) ! we cannot know the size ( it will be up to the user to 321 ! Have the correct size available 322 integer :: i,j,DA 323 324 DA = ubound(array,dim=1) 325 326 ! Everything else has to be checked... 327 UNIQ(1) = array(1) 328 j = 1 329 do i = 2 , DA 330 ! Apparently, inserting a do-loop is worse than count? 331 ! This makes me worried.... ??? 332! do ii = 1 , j 333! if ( UNIQ(ii) == array(i) ) cycle find_num 334! end do 335! j = j + 1 336! UNIQ(j) = array(i) 337 338 if ( count(UNIQ(1:j) == array(i)) == 0 ) then 339 j = j + 1 340 UNIQ(j) = array(i) 341 end if 342 343 end do 344 345 end function UNIQ 346 347! Returns an integer array sorted 348! It will contain dublicates if encountered. And hence sorting will amount to 349! arrays (/1,2,2,3,4,4,5/) for example. 350! This sorting routine has been optimized for consecutive segments 351! in the original array, hence, sorting on arrays with: 352! [1,2,3,4,25,26,27,28,15,16,17] are VERY fast! 353 pure function FSORT(array) result(SO) 354 integer, intent(in) :: array(:) 355 integer :: SO(ubound(array,dim=1)) 356 integer :: i,j,DA,h,FM 357 358 DA = ubound(array,dim=1) 359 360 ! The first entry is of course not found anywhere else 361 if ( DA == 0 ) then 362 SO = array 363 return 364 else if ( DA == 1 ) then 365 SO = array 366 return 367 end if 368 369 ! Everything else has to be checked... 370 ! This sorting routine is very efficient, and it works 371 ! Consider to change this function to a quick-sort algorithm... 372 SO(1) = array(1) 373 i = 2 374 sort_loop: do while ( i <= DA ) 375 if ( array(i) <= SO(1) ) then 376 !print '(a,1000(tr1,i0))','F1',SO(1:i-1),array(i) 377 ! put in front of the sorted array 378 call insert_front(DA,SO,array,i,FM) 379 i = i + FM 380 !print '(a,1000(tr1,i0))','F2',SO(1:i-1) 381 else if ( SO(i-1) <= array(i) ) then 382 !print '(a,1000(tr1,i0))','B1',SO(1:i-1),array(i) 383 call insert_back(DA,SO,array,i,FM) 384 i = i + FM 385 !print '(a,1000(tr1,i0))','B2',SO(1:i-1) 386 else if ( i-1 < 35 ) then 387 !print '(a,1000(tr1,i0))','M1',SO(1:i-1),array(i) 388 ! We assume that for array segments below 35 elements 389 ! it will be faster to do array search 390 expSearch: do j = 2 , i - 1 391 ! It will always be better to search for the end 392 ! of the overlapping region. i.e. less elements to move 393 if ( array(i) <= SO(j) ) then 394 h = j 395 call insert_mid(DA,SO,array,h,i,FM) 396 i = i + FM 397 exit expSearch ! exit the loop 398 end if 399 end do expSearch 400 !print '(a,1000(tr1,i0))','M2',SO(1:i-1) 401 402 else 403 !print '(a,1000(tr1,i0))','S1',SO(1:i-1),array(i) 404 405 ! search using SFIND, 406 ! We are taking advantage that both SO(1) and SO(i-1) 407 ! has been checked, hence the -1 408 j = SFIND(SO(2:i-1),array(i),NEAREST=+1) + 1 409 410 ! Insert directly, we have found what we searched for 411 call insert_mid(DA,SO,array,j,i,FM) 412 i = i + FM 413 414 !print '(a,1000(tr1,i0))','S2',SO(1:i-1) 415 end if 416 end do sort_loop 417 418 !do i = 1 , DA -1 419 ! if ( so(i) > so(i+1) )then 420 ! print *,i,DA,so(i),so(i+1) 421 ! print *,sum(so),sum(array) 422 ! call die('wrong sort') 423 ! end if 424 !end do 425 426 !if ( sum(so) /= sum(array) ) then 427 ! print *,sum(so),sum(array) 428 ! call die('wrong sort') 429 !end if 430 431 contains 432 433 pure subroutine insert_mid(DA,SO,array,sF,sA,P) 434 integer, intent(in) :: DA 435 integer, intent(in out) :: SO(DA) 436 integer, intent(in) :: array(DA) 437 ! The place where we found a SO(sF) <= array(sA) 438 integer, intent(in out) :: sF 439 integer, intent(in) :: sA ! The current reached iteration in array 440 integer, intent(out) :: P ! the number of inserted values 441 442 ! The last insertion point 443 integer :: lA 444 445 ! First we will skip to the last non-SAME value 446 ! I.e. where we can do the insert in SO 447 do while ( sF < sA - 1 .and. SO(sF) == array(sA) ) 448 sF = sF + 1 449 end do 450 451 ! Now SO(sF) < array(sA) 452 if ( sF >= sA ) then 453 call insert_back(DA,SO,array,sA,P) 454 return 455 end if 456 457!******* OLD 458 ! We wish to find the last element of array 459 ! we can move into the sort'ed array 460! lA = sA + 1 461 ! We know we are in the middle of the sort array 462 ! hence, we can exploit the next element in the sorted 463 ! array 464! do while ( SO(sF-1) <= array(lA-1) .and. & 465! array(lA) < SO(sF) .and. & 466! array(lA-1) <= array(lA) .and. & 467! lA <= DA ) ! must be in the range [sF;sA-1] 468! lA = lA + 1 469! end do 470! do while ( SO(sF) == array(lA) .and. lA <= DA ) 471 ! must be in the range [sF;sA-1] 472! lA = lA + 1 473! end do 474 475 476 lA = sA + 1 477 ! We know we are in the middle of the sort array 478 ! hence, we can exploit the next element in the sorted 479 ! array 480 do while ( lA <= DA ) 481 ! If the previous SO value is larger than the insertion 482 if ( SO(sF-1) > array(lA - 1) ) exit 483 ! If the array is not consecutive 484 if ( array(lA-1) > array(lA) ) exit 485 ! If the insertion point array is not consecutive 486 if ( array(lA) > SO(sF) ) exit 487 ! We need to ensure an overcount of 1 488 lA = lA + 1 489 end do 490 !if ( lA <= DA ) then 491 ! do while ( SO(sF) == array(lA) .and. lA <= DA ) 492 ! must be in the range [sF;sA-1] 493 ! lA = lA + 1 494 !end do 495 496 ! The number of elements we are pushing in 497 P = lA - sA 498 ! We have "overcounted" 499 !lA = lA - 1 500 501 !print '(6(tr1,a,tr1,i4))', & 502 ! 'SO |-1| =',SO(sF-1), & 503 ! 'SO |0| =',SO(sF), & 504 ! 'SO |+1| =',SO(sF+1) 505 !print '(6(tr1,a,tr1,i4))','P===',P,& 506 ! 'LH1',sA+P-1-(sF+P),& 507 ! 'RH1',sA-1-sF,& 508 ! 'LH2',sF+P-1-sF,& 509 ! 'RH2',sA+P-1-sA 510 !print '(6(tr1,i4))',array(sA:sA+5) 511 512 ! Copy the mid to the front of the SO 513 SO(sF+P:sA+P-1) = SO(sF:sA-1) 514 SO(sF:sF+P-1) = array(sA:sA+P-1) 515 516 end subroutine insert_mid 517 518 pure subroutine insert_front(DA,SO,array,sA,P) 519 integer, intent(in) :: DA 520 integer, intent(in out) :: SO(DA) 521 integer, intent(in) :: array(DA) 522 integer, intent(in) :: sA 523 integer, intent(out) :: P ! The number of Pasted values 524 525 ! The last insertion point 526 integer :: lA, i 527 528 i = sA + 1 529 do lA = i , DA 530 ! if the previous element is larger than the current element 531 if ( array(lA-1) >= array(lA) ) exit 532 ! if the checked point is larger than the insertion point 533 if ( array(lA) >= SO(1) ) exit 534 end do 535 i = lA 536 do lA = i , DA 537 if ( array(lA) /= SO(1) ) exit 538 end do 539 540 541! lA = sA + 1 542 ! Take all the values which are smaller than SO(1) 543! do while ( array(lA-1) < array(lA) .and. & 544! array(lA) < SO(1) .and. & 545! lA <= DA ) 546! lA = lA + 1 547! end do 548 ! Take all the values which are EQUAL to SO(1) 549! do while ( array(lA) == SO(1) .and. lA <= DA ) 550! lA = lA + 1 551! end do 552 553 ! Number of points found in the sort routine 554 P = lA - sA 555 ! Correct the overstepping (remark the above line counts correctly!) 556 !lA = lA - 1 557 !print '(6(tr1,a,tr1,i4))',& 558 ! 'LH1',P+sA-1-(1+P),& 559 ! 'RH1',sA-1-1,& 560 ! 'LH2',P-1,& 561 ! 'RH2',lA-sA 562 563 ! Copy over the values 564 SO(P+1:P+sA-1) = SO(1:sA-1) 565 SO(1:P) = array(sA:lA-1) 566 567 end subroutine insert_front 568 569 pure subroutine insert_back(DA,SO,array,sA,P) 570 integer, intent(in) :: DA 571 integer, intent(in out) :: SO(DA) 572 integer, intent(in) :: array(DA) 573 integer, intent(in) :: sA 574 integer, intent(out) :: P ! 575 576 ! The last insertion point 577 integer :: lA,i 578 579 lA = sA + 1 580 ! Step until SO(sA-1) /= array(lA) 581 do lA = sA + 1 , DA 582 if ( SO(sA-1) /= array(lA-1) ) exit 583 end do 584! OLD 585! do while ( SO(sA-1) == array(lA-1) .and. lA <= DA ) 586! lA = lA + 1 587! end do 588 589 ! Step until the last element of SO, SO(sA-1), is not 590 ! smaller than the array value 591 i = lA 592 do lA = i , DA 593 if ( SO(sA-1) >= array(lA-1) ) exit 594 if ( array(lA-1) >= array(lA) ) exit 595 end do 596! do while ( SO(sA-1) < array(lA-1) .and. & 597! array(lA-1) < array(lA) .and. & ! asserts the elements we choose are consecutive 598! lA <= DA ) 599! lA = lA + 1 600! end do 601 602 ! The number of elements we are pushing in 603 P = lA - sA 604 !print '(6(tr1,a,tr1,i4))',& 605 ! 'L/R H',sA+P-1-sA 606 SO(sA:sA+P-1) = array(sA:sA+P-1) 607 608 end subroutine insert_back 609 610 end function FSORT 611 612 ! Returns an integer array sorted 613! It will contain dublicates if encountered. And hence sorting will amount to 614! arrays (/1,2,2,3,4,4,5/) for example. 615! This sorting routine has been optimized for consecutive segments 616! in the original array, hence, sorting on arrays with: 617! [1,2,3,4,25,26,27,28,15,16,17] are VERY fast! 618 pure subroutine SORT(DA, array, so) 619 integer, intent(in) :: DA 620 integer, intent(in) :: array(DA) 621 integer, intent(out) :: SO(DA) 622 integer :: i, j, h, FM 623 624 ! No sorting 625 if ( DA == 0 ) return 626 627 SO(1) = array(1) 628 if ( DA == 1 ) return 629 630 ! Everything else has to be checked... 631 ! This sorting routine is very efficient, and it works 632 ! Consider to change this function to a quick-sort algorithm... 633 i = 2 634 sort_loop: do while ( i <= DA ) 635 if ( array(i) <= SO(1) ) then 636 !print '(a,1000(tr1,i0))','F1',SO(1:i-1),array(i) 637 ! put in front of the sorted array 638 call insert_front(DA,SO,array,i,FM) 639 i = i + FM 640 !print '(a,1000(tr1,i0))','F2',SO(1:i-1) 641 else if ( SO(i-1) <= array(i) ) then 642 !print '(a,1000(tr1,i0))','B1',SO(1:i-1),array(i) 643 call insert_back(DA,SO,array,i,FM) 644 i = i + FM 645 !print '(a,1000(tr1,i0))','B2',SO(1:i-1) 646 else if ( i-1 < 35 ) then 647 !print '(a,1000(tr1,i0))','M1',SO(1:i-1),array(i) 648 ! We assume that for array segments below 35 elements 649 ! it will be faster to do array search 650 expSearch: do j = 2 , i - 1 651 ! It will always be better to search for the end 652 ! of the overlapping region. i.e. less elements to move 653 if ( array(i) <= SO(j) ) then 654 h = j 655 call insert_mid(DA,SO,array,h,i,FM) 656 i = i + FM 657 exit expSearch ! exit the loop 658 end if 659 end do expSearch 660 !print '(a,1000(tr1,i0))','M2',SO(1:i-1) 661 662 else 663 !print '(a,1000(tr1,i0))','S1',SO(1:i-1),array(i) 664 665 ! search using SFIND, 666 ! We are taking advantage that both SO(1) and SO(i-1) 667 ! has been checked, hence the -1 668 j = SFIND(SO(2:i-1),array(i),NEAREST=+1) + 1 669 670 ! Insert directly, we have found what we searched for 671 call insert_mid(DA,SO,array,j,i,FM) 672 i = i + FM 673 674 !print '(a,1000(tr1,i0))','S2',SO(1:i-1) 675 end if 676 end do sort_loop 677 678 !do i = 1 , DA -1 679 ! if ( so(i) > so(i+1) )then 680 ! print *,i,DA,so(i),so(i+1) 681 ! print *,sum(so),sum(array) 682 ! call die('wrong sort') 683 ! end if 684 !end do 685 686 !if ( sum(so) /= sum(array) ) then 687 ! print *,sum(so),sum(array) 688 ! call die('wrong sort') 689 !end if 690 691 contains 692 693 pure subroutine insert_mid(DA,SO,array,sF,sA,P) 694 integer, intent(in) :: DA 695 integer, intent(in out) :: SO(DA) 696 integer, intent(in) :: array(DA) 697 ! The place where we found a SO(sF) <= array(sA) 698 integer, intent(in out) :: sF 699 integer, intent(in) :: sA ! The current reached iteration in array 700 integer, intent(out) :: P ! the number of inserted values 701 702 ! The last insertion point 703 integer :: lA, i 704 705 ! First we will skip to the last non-SAME value 706 ! I.e. where we can do the insert in SO 707 do while ( sF < sA - 1 .and. SO(sF) == array(sA) ) 708 sF = sF + 1 709 end do 710 711 ! Now SO(sF) < array(sA) 712 if ( sF >= sA ) then 713 call insert_back(DA,SO,array,sA,P) 714 return 715 end if 716 717!******* OLD 718 ! We wish to find the last element of array 719 ! we can move into the sort'ed array 720! lA = sA + 1 721 ! We know we are in the middle of the sort array 722 ! hence, we can exploit the next element in the sorted 723 ! array 724! do while ( SO(sF-1) <= array(lA-1) .and. & 725! array(lA) < SO(sF) .and. & 726! array(lA-1) <= array(lA) .and. & 727! lA <= DA ) ! must be in the range [sF;sA-1] 728! lA = lA + 1 729! end do 730! do while ( SO(sF) == array(lA) .and. lA <= DA ) 731 ! must be in the range [sF;sA-1] 732! lA = lA + 1 733! end do 734 735 736 lA = sA + 1 737 ! We know we are in the middle of the sort array 738 ! hence, we can exploit the next element in the sorted 739 ! array 740 do while ( lA <= DA ) 741 ! If the previous SO value is larger than the insertion 742 if ( SO(sF-1) > array(lA - 1) ) exit 743 ! If the array is not consecutive 744 if ( array(lA-1) > array(lA) ) exit 745 ! If the insertion point array is not consecutive 746 if ( array(lA) > SO(sF) ) exit 747 ! We need to ensure an overcount of 1 748 lA = lA + 1 749 end do 750 !if ( lA <= DA ) then 751 ! do while ( SO(sF) == array(lA) .and. lA <= DA ) 752 ! must be in the range [sF;sA-1] 753 ! lA = lA + 1 754 !end do 755 756 ! The number of elements we are pushing in 757 P = lA - sA 758 ! We have "overcounted" 759 !lA = lA - 1 760 761 !print '(6(tr1,a,tr1,i4))', & 762 ! 'SO |-1| =',SO(sF-1), & 763 ! 'SO |0| =',SO(sF), & 764 ! 'SO |+1| =',SO(sF+1) 765 !print '(6(tr1,a,tr1,i4))','P===',P,& 766 ! 'LH1',sA+P-1-(sF+P),& 767 ! 'RH1',sA-1-sF,& 768 ! 'LH2',sF+P-1-sF,& 769 ! 'RH2',sA+P-1-sA 770 !print '(6(tr1,i4))',array(sA:sA+5) 771 772 ! Copy the mid to the front of the SO 773 if ( sF + P <= sA - 1 ) then 774 SO(sF+P:sA+P-1) = SO(sF:sA-1) 775 else 776 do i = sF , sA - 1 777 SO(i+P) = SO(i) 778 end do 779 end if 780 !SO(sF:sF+P-1) = array(sA:sA+P-1) 781 do i = 0 , P - 1 782 SO(sF+i) = array(sA+i) 783 end do 784 785 end subroutine insert_mid 786 787 pure subroutine insert_front(DA,SO,array,sA,P) 788 integer, intent(in) :: DA 789 integer, intent(in out) :: SO(DA) 790 integer, intent(in) :: array(DA) 791 integer, intent(in) :: sA 792 integer, intent(out) :: P ! The number of Pasted values 793 794 ! The last insertion point 795 integer :: lA, i 796 797 i = sA + 1 798 do lA = i , DA 799 ! if the previous element is larger than the current element 800 if ( array(lA-1) >= array(lA) ) exit 801 ! if the checked point is larger than the insertion point 802 if ( array(lA) >= SO(1) ) exit 803 end do 804 i = lA 805 do lA = i , DA 806 if ( array(lA) /= SO(1) ) exit 807 end do 808 809 810! lA = sA + 1 811 ! Take all the values which are smaller than SO(1) 812! do while ( array(lA-1) < array(lA) .and. & 813! array(lA) < SO(1) .and. & 814! lA <= DA ) 815! lA = lA + 1 816! end do 817 ! Take all the values which are EQUAL to SO(1) 818! do while ( array(lA) == SO(1) .and. lA <= DA ) 819! lA = lA + 1 820! end do 821 822 ! Number of points found in the sort routine 823 P = lA - sA 824 ! Correct the overstepping (remark the above line counts correctly!) 825 !lA = lA - 1 826 !print '(6(tr1,a,tr1,i4))',& 827 ! 'LH1',P+sA-1-(1+P),& 828 ! 'RH1',sA-1-1,& 829 ! 'LH2',P-1,& 830 ! 'RH2',lA-sA 831 832 ! Copy over the values 833 if ( P + 1 <= sA - 1 ) then 834 SO(P+1:P+sA-1) = SO(1:sA-1) 835 else 836 do i = 1 , sA - 1 837 SO(P+i) = SO(i) 838 end do 839 end if 840 !SO(1:P) = array(sA:lA-1) 841 do i = 1 , P 842 SO(i) = array(sA+i-1) 843 end do 844 845 end subroutine insert_front 846 847 pure subroutine insert_back(DA,SO,array,sA,P) 848 integer, intent(in) :: DA 849 integer, intent(in out) :: SO(DA) 850 integer, intent(in) :: array(DA) 851 integer, intent(in) :: sA 852 integer, intent(out) :: P ! 853 854 ! The last insertion point 855 integer :: lA,i 856 857 lA = sA + 1 858 ! Step until SO(sA-1) /= array(lA) 859 do lA = sA + 1 , DA 860 if ( SO(sA-1) /= array(lA-1) ) exit 861 end do 862! OLD 863! do while ( SO(sA-1) == array(lA-1) .and. lA <= DA ) 864! lA = lA + 1 865! end do 866 867 ! Step until the last element of SO, SO(sA-1), is not 868 ! smaller than the array value 869 i = lA 870 do lA = i , DA 871 if ( SO(sA-1) >= array(lA-1) ) exit 872 if ( array(lA-1) >= array(lA) ) exit 873 end do 874! do while ( SO(sA-1) < array(lA-1) .and. & 875! array(lA-1) < array(lA) .and. & ! asserts the elements we choose are consecutive 876! lA <= DA ) 877! lA = lA + 1 878! end do 879 880 ! The number of elements we are pushing in 881 P = lA - sA 882 !print '(6(tr1,a,tr1,i4))',& 883 ! 'L/R H',sA+P-1-sA 884 do i = 0 , P - 1 885 SO(sA+i) = array(sA+i) 886 end do 887 888 end subroutine insert_back 889 890 end subroutine SORT 891 892 subroutine index_sort_heap_i4( N, A, IDX ) 893! ******************************************************************* 894! SUBROUTINE index_sort_heap_i4( N, A, IDX ) 895! Makes an index table of increasing array A, with size N 896! using the heapsort algorithm. 897! Written by J.M.Soler, May.2015 898! Re-written by N.R.Papior, Mar.2018, taken from ORDIX (sorting.f) 899! *************** INPUT ********************************************* 900! INTEGER N : Dimensions of array A 901! integer A(N) : Array with the values to be ordered 902! *************** OUTPUT ******************************************** 903! INTEGER IDX(N) : Array which gives the increasing order of A(I): 904! A(IDX(I)) .LE. A(IDX(I+1)) 905! *************** ALGORITHM ***************************************** 906! A hierarchical 'family tree' (heap) is generated, with each parent 907! k older than its two children (2*k and 2*k+1). Persons k with k>np 908! are children (np = highest power of 2 smaller than n). Persons with 909! np/2<k<=np are parents (but only those with k<=n/2 have actually 910! one or two children). Persons np/4<k<=np/2 are grandparents, etc. 911! Then, the person with k=1 (the oldest) is removed and the tree is 912! reconstructed. This is iterated until all members are picked. 913! Ref: W.H.Press et al. Numerical Recipes, Cambridge Univ. Press. 914! ******************************************************************* 915 916 integer, intent(in) :: N ! Dimensions of array x 917 integer, intent(in) :: A(n) ! Array with the values to be ordered 918 ! will order against the first element (m) 919 integer, intent(inout) :: idx(n) ! Increasing order of A(:) 920 921 integer:: k, nFamily, parent 922 923 ! Construct the heap (family tree) 924 idx = (/(k,k=1,n)/) ! initial array order, to be modified 925 926 927 ! Swap to create the actual parent tree 928 nFamily = n ! number of persons in the family tree 929 do parent = n/2,1,-1 ! sift 'parents' down the tree 930 call siftDown(parent) ! siftDown inherits age and indx arrays 931 end do 932 933 ! Reduce the tree size, retiring its succesive patriarchs (first element) 934 do nFamily = n-1,1,-1 ! nFamily is the new size of the tree 935 call swap( idx(1), idx(nFamily+1) ) ! swap patriarch and youngest child 936 call siftDown(1) ! now recolocate child in tree 937 end do 938 939 contains 940 941 subroutine siftDown( person ) ! place person in family tree 942 integer, intent(in) :: person 943 944 ! Inherited from ordix: age(:), ageTol, idx(:), nFamily 945 integer:: child, sw, parent 946 947 ! initialize 948 parent = person ! assume person is a parent 949 child = 2 * parent ! first child of parent 950 sw = parent ! sorted swap child 951 do while ( child <= nFamily ) ! iterate the sift-down process 952 ! check current child 953 if ( A(idx(sw)) < A(idx(child)) ) then 954 sw = child 955 end if 956 ! Check neighbouring child 957 if ( child < nFamily ) then 958 if ( A(idx(sw)) < A(idx(child+1)) ) then 959 sw = child + 1 960 end if 961 end if 962 if ( sw == parent ) then 963 exit ! break 964 else 965 call swap( idx(parent) , idx(sw) ) 966 ! update for next sift 967 parent = sw 968 child = sw * 2 969 end if 970 end do 971 972 end subroutine siftDown 973 974 subroutine swap(i,j) ! exchange integers i and j 975 integer:: i,j,k 976 k = i 977 i = j 978 j = k 979 end subroutine swap 980 981 end subroutine index_sort_heap_i4 982 983 984 recursive pure subroutine sort_quick(n, array) 985 integer, intent(in) :: n 986 integer, intent(inout) :: array(n) 987 988 integer :: div 989 990 if ( n <= 1 ) return 991 992 ! Retrieve the partition ID 993 call partition(n, array, div) 994 call sort_quick(div-1, array(1)) 995 call sort_quick(n-div+1, array(div)) 996 997 contains 998 999 ! Partition an array by swapping 1000 pure subroutine partition(n, array, mark) 1001 integer, intent(in) :: n 1002 integer, intent(inout) :: array(n) 1003 integer, intent(out) :: mark 1004 1005 integer :: i, j, x, t 1006 1007 ! Simple check for n == 2 1008 ! Note that n == 1 will NEVER happen because 1009 ! of the sort_quick check 1010 if ( n == 2 ) then 1011 if ( array(2) < array(1) ) then 1012 t = array(2) 1013 array(2) = array(1) 1014 array(1) = t 1015 end if 1016 mark = 2 1017 return 1018 end if 1019 1020 i = n / 2 1021 ! Find the median of the searched elements 1022 if ( array(1) < array(n) ) then 1023 if ( array(n) < array(i) ) then 1024 x = array(n) 1025 else if ( array(1) < array(i) ) then 1026 x = array(i) 1027 else 1028 x = array(1) 1029 end if 1030 else !if ( array(n) < array(1) ) 1031 if ( array(1) < array(i) ) then 1032 x = array(1) 1033 else if ( array(n) < array(i) ) then 1034 x = array(i) 1035 else 1036 x = array(n) 1037 end if 1038 end if 1039 1040 ! Starting counters.. 1041 i = 0 1042 j = n + 1 1043 1044 do 1045 ! find point from below which is 1046 ! above median 1047 j = j - 1 1048 do while ( j > 0 ) 1049 if ( array(j) <= x ) exit 1050 j = j - 1 1051 end do 1052 1053 i = i + 1 1054 do while ( i < n ) 1055 if ( array(i) >= x ) exit 1056 i = i + 1 1057 end do 1058 1059 if ( i < j ) then 1060 ! exchange array(i) and array(j) 1061 t = array(i) 1062 array(i) = array(j) 1063 array(j) = t 1064 else if ( i == j ) then 1065 mark = i + 1 1066 return 1067 else 1068 mark = i 1069 return 1070 end if 1071 end do 1072 1073 end subroutine partition 1074 1075 end subroutine sort_quick 1076 1077! ***************** old CORRECT version ****************** 1078! Returns an integer array sorted 1079! It will contain dublicates if encountered. And hence sorting will amount to 1080! arrays (/1,2,2,3,4,4,5/) for example. 1081! function SORT(array) 1082! integer, intent(in) :: array(:) 1083! integer :: SORT(ubound(array,dim=1)) 1084! integer :: i,j,sa 1085! 1086! sa = size(array) 1087! 1088! ! The first entry is of course not found anywhere else 1089! if ( sa == 0 ) then 1090! SORT = array 1091! return 1092! else if ( sa == 1 ) then 1093! SORT = array 1094! return 1095! end if 1096! SORT(1) = array(1) 1097! sort_loop: do i = 2 , sa 1098! do j = 1, i-1 1099! It will always be better to search for the end 1100! of the overlapping region. i.e. less elements to move 1101! if ( array(i) < SORT(j) ) then 1102! SORT(j+1:i) = SORT(j:i-1) 1103! SORT(j) = array(i) 1104! cycle sort_loop 1105! end if 1106! end do 1107! SORT(i) = array(i) 1108! end do sort_loop 1109! end function SORT 1110 1111 1112! Functions for returning the Euclediean norm of a vector 1113! This takes two function interfaces: 1114! VNORM(vec(:,:)) which will compute the length of each column and return 1115! a vector having length ubound(vec(:,:),dim=2), i.e. 1116! vector length of fastests index 1117 pure function VNORM_sp_1(vec) result(norm) 1118 real(sp), intent(in) :: vec(:) 1119 real(sp) :: norm 1120 integer :: i 1121 norm = vec(1) * vec(1) 1122 do i = 2 , ubound(vec,dim=1) 1123 norm = norm + vec(i) * vec(i) 1124 end do 1125 norm = sqrt(norm) 1126 end function VNORM_sp_1 1127 pure function VNORM_sp_2(vec) result(norm) 1128 real(sp), intent(in) :: vec(:,:) 1129 real(sp) :: norm(ubound(vec,dim=2)) 1130 integer :: i 1131 do i = 1 , ubound(vec,dim=2) 1132 norm(i) = VNORM(vec(:,i)) 1133 end do 1134 end function VNORM_sp_2 1135 pure function VNORM_dp_1(vec) result(norm) 1136 real(dp), intent(in) :: vec(:) 1137 real(dp) :: norm 1138 integer :: i 1139 norm = vec(1) * vec(1) 1140 do i = 2 , ubound(vec,dim=1) 1141 norm = norm + vec(i) * vec(i) 1142 end do 1143 norm = dsqrt(norm) 1144 end function VNORM_dp_1 1145 pure function VNORM_dp_2(vec) result(norm) 1146 real(dp), intent(in) :: vec(:,:) 1147 real(dp) :: norm(ubound(vec,dim=2)) 1148 integer :: i 1149 do i = 1 , ubound(vec,dim=2) 1150 norm(i) = VNORM(vec(:,i)) 1151 end do 1152 end function VNORM_dp_2 1153 pure function VNORM_cp_1(vec) result(norm) 1154 complex(sp), intent(in) :: vec(:) 1155 real(sp) :: norm 1156 integer :: i 1157 norm = conjg(vec(1)) * vec(1) 1158 do i = 2 , ubound(vec,dim=1) 1159 norm = norm + conjg(vec(i)) * vec(i) 1160 end do 1161 norm = sqrt(norm) 1162 end function VNORM_cp_1 1163 pure function VNORM_cp_2(vec) result(norm) 1164 complex(sp), intent(in) :: vec(:,:) 1165 real(sp) :: norm(ubound(vec,dim=2)) 1166 integer :: i 1167 do i = 1 , ubound(vec,dim=2) 1168 norm(i) = VNORM(vec(:,i)) 1169 end do 1170 end function VNORM_cp_2 1171 pure function VNORM_zp_1(vec) result(norm) 1172 complex(dp), intent(in) :: vec(:) 1173 real(dp) :: norm 1174 integer :: i 1175 norm = dconjg(vec(1)) * vec(1) 1176 do i = 2 , ubound(vec,dim=1) 1177 norm = norm + dconjg(vec(i)) * vec(i) 1178 end do 1179 norm = dsqrt(norm) 1180 end function VNORM_zp_1 1181 pure function VNORM_zp_2(vec) result(norm) 1182 complex(dp), intent(in) :: vec(:,:) 1183 real(dp) :: norm(ubound(vec,dim=2)) 1184 integer :: i 1185 do i = 1 , ubound(vec,dim=2) 1186 norm(i) = VNORM(vec(:,i)) 1187 end do 1188 end function VNORM_zp_2 1189 1190 pure function ANGLE_sp(v1,v2) result(ang) 1191 real(sp), intent(in) :: v1(:), v2(:) 1192 real(sp) :: ang 1193 ang = acos( sum(v1 * v2) / (VNORM(v1)+VNORM(v2)) ) 1194 end function ANGLE_sp 1195 pure function ANGLE_dp(v1,v2) result(ang) 1196 real(dp), intent(in) :: v1(:), v2(:) 1197 real(dp) :: ang 1198 ang = acos( sum(v1 * v2) / (VNORM(v1)+VNORM(v2)) ) 1199 end function ANGLE_dp 1200 1201 1202! We add an algorithm for search for SORTED entries 1203! This will relieve a lot of programming in other segments 1204! of the codes. 1205! What it does is the following: 1206! 1. Get a sorted array of some size 1207! 2. Get a value which should be found in the array 1208! 3. Returns the index of the value found in the array 1209! 4. If the value is not found, return 0 1210! 5. If something went wrong, return -1 1211 pure function SORTED_FIND(array,val) result(idx) 1212 integer, intent(in) :: array(:) 1213 integer, intent(in) :: val 1214 ! Used internal variables 1215 integer :: DA,h,FM, idx 1216 1217 ! Retrieve the size of the array we search in 1218 DA = ubound(array,dim=1) 1219 1220 ! Initialize to default value 1221 idx = 0 1222 if ( DA == 0 ) return 1223 1224 ! The two easiest cases, i.e. they are not in the array... 1225 if ( val < array(1) ) then 1226 return 1227 else if ( val == array(1) ) then 1228 idx = 1 1229 return 1230 else if ( array(DA) < val ) then 1231 return 1232 else if ( val == array(DA) ) then 1233 idx = DA 1234 return 1235 end if 1236 1237 ! If it is size 1 or 2, we can immediately return, 1238 ! we have already checked the first/last index 1239 if ( DA <= 2 ) return 1240 1241 ! An *advanced* search algorithm... 1242 1243 ! Search the sorted array for an entry 1244 ! We know it *must* have one 1245 h = DA / 2 1246 idx = h ! Start in the middle 1247 ! The integer correction (due to round of errors when 1248 ! calculating the new half... 1249 FM = MOD(h,2) 1250 do while ( h > 0 ) ! While we are still searching... 1251 1252 if ( h >= 3 ) then 1253 h = h + FM 1254 ! This will only add 1 every other time 1255 FM = MOD(h,2) 1256 end if 1257 1258 ! integer division is fast 1259 h = h / 2 1260 1261 if ( val < array(idx) ) then 1262 ! the value we search for is smaller than 1263 ! the current checked value, hence we step back 1264 !print *,'stepping down',i,h 1265 idx = idx - h 1266 else if ( array(idx) < val ) then 1267 ! the value we search for is larger than 1268 ! the current checked value, hence we step forward 1269 !print *,'stepping up',i,h 1270 idx = idx + h 1271 else 1272 !print *,'found',i 1273 ! We know EXACTLY where we are... 1274 return 1275 ! We found it!!! 1276 end if 1277 end do 1278 1279 ! We need to ensure a range to search in 1280 ! The missing integers are *only* necesseary when 1281 ! the search pattern is in the same direction. 1282 ! This can easily be verified... 1283 h = 1 + FM 1284 1285 ! The missing integer count ensures the correct range 1286 FM = max(idx - h, 1 ) 1287 h = min(idx + h, DA) 1288 1289 ! The index will *most* likely be close to 'i' 1290 ! Hence we start by searching around it 1291 ! However, we know that val /= array(i) 1292 1293 do idx = FM, h 1294 if ( val == array(idx) ) return 1295 end do 1296 1297 ! Default value is *not found* 1298 idx = 0 1299 1300 end function SORTED_FIND 1301 1302 pure function SORTED_BINARY_FIND(array,val) result(idx) 1303 integer, intent(in) :: array(:) 1304 integer, intent(in) :: val 1305 1306 ! Used internal variables 1307 integer :: idx 1308 integer :: small, large, d 1309 1310 ! Retrieve the size of the array we search in 1311 d = ubound(array,dim=1) 1312 1313 ! Initialize to default value 1314 idx = 0 1315 if ( d == 0 ) return 1316 1317 ! The two easiest cases, i.e. they are not in the array... 1318 if ( val < array(1) ) then 1319 return 1320 else if ( val == array(1) ) then 1321 idx = 1 1322 return 1323 else if ( array(d) < val ) then 1324 return 1325 else if ( val == array(d) ) then 1326 idx = d 1327 return 1328 end if 1329 1330 ! If it is size 1 or 2, we can immediately return, 1331 ! we have already checked the first/last index 1332 if ( d <= 2 ) return 1333 1334 small = 1 1335 large = d 1336 do while ( small + 1 < large ) 1337 idx = (small + large) / 2 1338 if ( array(idx) < val ) then 1339 small = idx 1340 else 1341 large = idx 1342 end if 1343 end do 1344 if ( array(idx) == val ) return 1345 do idx = small, large 1346 if ( array(idx) == val ) return 1347 end do 1348 idx = 0 1349 1350 end function SORTED_BINARY_FIND 1351 1352 pure recursive function SORTED_REC_FIND(array,val) result(idx) 1353 integer, intent(in) :: array(:) 1354 integer, intent(in) :: val 1355 ! Used internal variables 1356 integer :: DA, idx, tmp 1357 1358 ! Retrieve the size of the array we search in 1359 DA = ubound(array,dim=1) 1360 1361 ! Initialize to default value 1362 idx = 0 1363 if ( DA == 0 ) return 1364 1365 ! The two easiest cases, i.e. they are not in the array... 1366 if ( val < array(1) ) then 1367 return 1368 else if ( val == array(1) ) then 1369 idx = 1 1370 return 1371 else if ( array(DA) < val ) then 1372 return 1373 else if ( val == array(DA) ) then 1374 idx = DA 1375 return 1376 end if 1377 1378 ! If it is size 1, we can immediately return, 1379 ! we have already checked the first index 1380 if ( DA <= 2 ) then 1381 idx = 0 1382 return 1383 end if 1384 1385 idx = DA / 2 1386 if ( val <= array(idx) ) then 1387 idx = SORTED_REC_FIND(array(1:idx),val) 1388 else 1389 tmp = SORTED_REC_FIND(array(idx+1:),val) 1390 if ( tmp == 0 ) then 1391 idx = 0 1392 else 1393 idx = idx + tmp 1394 end if 1395 end if 1396 1397 end function SORTED_REC_FIND 1398 1399! We add an algorithm for search for SORTED entries 1400! This will relieve a lot of programming in other segments 1401! of the codes. 1402! What it does is the following: 1403! 1. Get a sorted array of some size 1404! 2. Get a value which should be found in the array 1405! 3. Returns the index of the value found in the array 1406! 4. If the value is not found, return 0 1407! 5. If something went wrong, return -1 1408 pure function SORTED_CLOSE(array,val,NEAREST) result(idx) 1409 integer, intent(in) :: array(:) 1410 integer, intent(in) :: val 1411 integer, intent(in) :: NEAREST 1412 ! Used internal variables 1413 integer :: DA,h,FM, idx 1414 1415 ! Retrieve the size of the array we search in 1416 DA = ubound(array,dim=1) 1417 1418 ! Initialize to default value 1419 idx = 0 1420 if ( DA == 0 ) return 1421 1422 ! The two easiest cases, i.e. they are not in the array... 1423 if ( val < array(1) ) then 1424 ! We need only handle if the user requests 1425 ! *closest above* 1426 if ( NEAREST > 0 ) idx = 1 1427 return 1428 else if ( val == array(1) ) then 1429 idx = 1 1430 return 1431 else if ( array(DA) < val ) then 1432 ! We need only handle if the user requests 1433 ! *closest below* 1434 if ( NEAREST < 0 ) idx = DA 1435 return 1436 else if ( val == array(DA) ) then 1437 idx = DA 1438 return 1439 end if 1440 1441 ! If it is size 1, we can immediately return, 1442 ! we have already checked the first index 1443 if ( DA == 1 ) return 1444 1445 ! An *advanced* search algorithm... 1446 1447 ! Search the sorted array for an entry 1448 ! We know it *must* have one 1449 h = DA / 2 1450 if ( h < 1 ) h = DA ! This ensures a correct handling 1451 idx = h ! Start in the middle 1452 ! The integer correction (due to round of errors when 1453 ! calculating the new half... 1454 FM = MOD(h,2) 1455 do while ( h > 1 ) ! While we are still searching... 1456 1457 if ( h >= 3 ) then 1458 h = h + FM 1459 ! This will only add 1 every other time 1460 FM = MOD(h,2) 1461 end if 1462 ! integer division is faster. :) 1463 h = h / 2 1464 1465 ! We must ensure to never step UNDER zero 1466 ! This we can only do by using "floor" 1467 !h = int(h*0.5_dp) 1468 1469 ! This makes sure that we 'track' the potential 1470 ! missing integer while performing 'floor' 1471 ! This missing integer, will only arise 1472 ! when the number is non-divisable by two 1473 !FM = FM + MOD(h,2) 1474 1475 if ( val < array(idx) ) then 1476 ! the value we search for is smaller than 1477 ! the current checked value, hence we step back 1478 !print *,'stepping down',i,h 1479 idx = idx - h 1480 else if ( array(idx) < val ) then 1481 ! the value we search for is larger than 1482 ! the current checked value, hence we step forward 1483 !print *,'stepping up',i,h 1484 idx = idx + h 1485 else 1486 !print *,'found',i 1487 ! We know EXACTLY where we are... 1488 return 1489 ! We found it!!! 1490 end if 1491 end do 1492 1493 ! We need to ensure a range to search in 1494 ! The missing integers are *only* necesseary when 1495 ! the search pattern is in the same direction. 1496 ! This can easily be verified... 1497 h = max(h,1) + FM + 1 1498 1499 ! The missing integer count ensures the correct range 1500 FM = max(idx - h, 1 ) 1501 h = min(idx + h, DA) 1502 1503 ! The index will *most* likely be close to 'i' 1504 ! Hence we start by searching around it 1505 ! However, we know that val /= array(i) 1506 1507 ! We need to determine the method of indexing 1508 if ( NEAREST < 0 ) then 1509 do idx = FM, h - 1 1510 if ( val == array(idx) ) then 1511 return 1512 else if ( val < array(idx+1) ) then 1513 return 1514 end if 1515 end do 1516 ! This checks for array(h) 1517 if ( val == array(h) ) then 1518 idx = h 1519 return 1520 end if 1521 else if ( NEAREST == 0 ) then 1522 do idx = FM, h 1523 if ( val == array(idx) ) return 1524 end do 1525 else if ( NEAREST > 0 ) then 1526 do idx = FM, h 1527 if ( val == array(idx) ) then 1528 return 1529 else if ( val < array(idx) .and. idx > 1 ) then 1530 if ( val > array(idx-1) ) return 1531 end if 1532 end do 1533 else 1534 ! ERROR 1535 idx = -1 1536 end if 1537 1538 ! Default value is *not found* 1539 idx = 0 1540 1541 end function SORTED_CLOSE 1542 1543 1544 ! Matrix operations missing 1545 ! The reason for choosing a subroutine for these 1546 ! are the direct impact on memory for very large matrices. 1547 ! This ensures direct writes, instead of temporary eye-arrays... 1548 subroutine EYE_i_2D(size,array,I) 1549 integer, intent(in) :: size 1550 integer, intent(out) :: array(size,size) 1551 integer, intent(in), optional :: I 1552 integer :: j, k, lI 1553 lI = 1 1554 if ( present(I) ) lI = I 1555!$OMP parallel do default(shared), private(k,j) 1556 do k = 1 , size 1557 do j = 1 , size 1558 array(j,k) = 0 1559 end do 1560 array(k,k) = lI 1561 end do 1562!$OMP end parallel do 1563 end subroutine EYE_i_2D 1564 subroutine EYE_sp_2D(size,array,I) 1565 integer, intent(in) :: size 1566 real(sp), intent(out) :: array(size,size) 1567 real(sp), intent(in), optional :: I 1568 real(sp) :: lI 1569 integer :: j, k 1570 lI = 1._sp 1571 if ( present(I) ) lI = I 1572!$OMP parallel do default(shared), private(k,j) 1573 do k = 1 , size 1574 do j = 1 , size 1575 array(j,k) = 0._sp 1576 end do 1577 array(k,k) = lI 1578 end do 1579!$OMP end parallel do 1580 end subroutine EYE_sp_2D 1581 subroutine EYE_dp_2D(size,array,I) 1582 integer, intent(in) :: size 1583 real(dp), intent(out) :: array(size,size) 1584 real(dp), intent(in), optional :: I 1585 real(dp) :: lI 1586 integer :: j, k 1587 lI = 1._dp 1588 if ( present(I) ) lI = I 1589!$OMP parallel do default(shared), private(k,j) 1590 do k = 1 , size 1591 do j = 1 , size 1592 array(j,k) = 0._dp 1593 end do 1594 array(k,k) = lI 1595 end do 1596!$OMP end parallel do 1597 end subroutine EYE_dp_2D 1598 subroutine EYE_cp_2D(size,array,I) 1599 integer, intent(in) :: size 1600 complex(sp), intent(out) :: array(size,size) 1601 complex(sp), intent(in), optional :: I 1602 complex(sp) :: lI 1603 integer :: j, k 1604 lI = cmplx(1._sp,0._sp) 1605 if ( present(I) ) lI = I 1606!$OMP parallel do default(shared), private(k,j) 1607 do k = 1 , size 1608 do j = 1 , size 1609 array(j,k) = 0._sp 1610 end do 1611 array(k,k) = lI 1612 end do 1613!$OMP end parallel do 1614 end subroutine EYE_cp_2D 1615 subroutine EYE_zp_2D(size,array,I) 1616 integer, intent(in) :: size 1617 complex(dp), intent(out) :: array(size,size) 1618 complex(dp), intent(in), optional :: I 1619 complex(dp) :: lI 1620 integer :: j, k 1621 lI = dcmplx(1._dp,0._dp) 1622 if ( present(I) ) lI = I 1623!$OMP parallel do default(shared), private(k,j) 1624 do k = 1 , size 1625 do j = 1 , size 1626 array(j,k) = dcmplx(0._dp,0._dp) 1627 end do 1628 array(k,k) = lI 1629 end do 1630!$OMP end parallel do 1631 end subroutine EYE_zp_2D 1632 1633 subroutine EYE_i_1D(size,array) 1634 integer, intent(in) :: size 1635 integer, intent(out) :: array(size*size) 1636 call EYE_i_2D(size,array) 1637 end subroutine EYE_i_1D 1638 subroutine EYE_sp_1D(size,array) 1639 integer, intent(in) :: size 1640 real(sp), intent(out) :: array(size*size) 1641 call EYE_sp_2D(size,array) 1642 end subroutine EYE_sp_1D 1643 subroutine EYE_dp_1D(size,array) 1644 integer, intent(in) :: size 1645 real(dp), intent(out) :: array(size*size) 1646 call EYE_dp_2D(size,array) 1647 end subroutine EYE_dp_1D 1648 subroutine EYE_cp_1D(size,array) 1649 integer, intent(in) :: size 1650 complex(sp), intent(out) :: array(size*size) 1651 call EYE_cp_2D(size,array) 1652 end subroutine EYE_cp_1D 1653 subroutine EYE_zp_1D(size,array) 1654 integer, intent(in) :: size 1655 complex(dp), intent(out) :: array(size*size) 1656 call EYE_zp_2D(size,array) 1657 end subroutine EYE_zp_1D 1658 1659 1660 subroutine TRANSPOSE_dp_2D(size,array) 1661 integer, intent(in) :: size 1662 real(dp), intent(inout) :: array(size,size) 1663 integer :: i, j 1664 real(dp) :: d 1665!$OMP parallel do default(shared), private(i,j,d) 1666 do i = 1 , size - 1 1667 do j = i + 1 , size 1668 d = array(j,i) 1669 array(j,i) = array(i,j) 1670 array(i,j) = d 1671 end do 1672 end do 1673!$OMP end parallel do 1674 end subroutine TRANSPOSE_dp_2D 1675 subroutine TRANSPOSE_zp_2D(size,array) 1676 integer, intent(in) :: size 1677 complex(dp), intent(inout) :: array(size,size) 1678 integer :: i, j 1679 complex(dp) :: z 1680!$OMP parallel do default(shared), private(i,j,z) 1681 do i = 1 , size - 1 1682 do j = i + 1 , size 1683 z = array(j,i) 1684 array(j,i) = array(i,j) 1685 array(i,j) = z 1686 end do 1687 end do 1688!$OMP end parallel do 1689 end subroutine TRANSPOSE_zp_2D 1690 1691 subroutine TRANSPOSE_dp_1D(size,array) 1692 integer, intent(in) :: size 1693 real(dp), intent(inout) :: array(size*size) 1694 call TRANSPOSE_dp_2D(size,array) 1695 end subroutine TRANSPOSE_dp_1D 1696 subroutine TRANSPOSE_zp_1D(size,array) 1697 integer, intent(in) :: size 1698 complex(dp), intent(inout) :: array(size*size) 1699 call TRANSPOSE_zp_2D(size,array) 1700 end subroutine TRANSPOSE_zp_1D 1701 1702 1703 function TRACE_dp_2D(size,array) result(T) 1704 integer, intent(in) :: size 1705 real(dp), intent(in) :: array(size,size) 1706 real(dp) :: T 1707 integer :: i 1708 T = 0._dp 1709!$OMP parallel do default(shared), private(i), reduction(+:T) 1710 do i = 1 , size 1711 T = T + array(i,i) 1712 end do 1713!$OMP end parallel do 1714 end function TRACE_dp_2D 1715 function TRACE_zp_2D(size,array) result(T) 1716 integer, intent(in) :: size 1717 complex(dp), intent(in) :: array(size,size) 1718 complex(dp) :: T 1719 integer :: i 1720 T = dcmplx(0._dp,0._dp) 1721!$OMP parallel do default(shared), private(i), reduction(+:T) 1722 do i = 1 , size 1723 T = T + array(i,i) 1724 end do 1725!$OMP end parallel do 1726 end function TRACE_zp_2D 1727 1728 function TRACE_dp_1D(size,array) result(T) 1729 integer, intent(in) :: size 1730 real(dp), intent(in) :: array(size*size) 1731 real(dp) :: T 1732 T = TRACE_dp_2D(size,array) 1733 end function TRACE_dp_1D 1734 function TRACE_zp_1D(size,array) result(T) 1735 integer, intent(in) :: size 1736 complex(dp), intent(in) :: array(size*size) 1737 complex(dp) :: T 1738 T = TRACE_zp_2D(size,array) 1739 end function TRACE_zp_1D 1740 1741 1742 ! Projections from one direction onto ND space 1743 pure function SPC_PROJ_sp(space,vin) result(vout) 1744 real(sp), intent(in) :: space(:,:), vin(:) 1745 real(sp) :: vout(size(vin)), tmp(size(vin)) 1746 integer :: i 1747 do i = 1 , size(vin) 1748 tmp(:) = space(:,i) / VNORM(space(:,i)) 1749 vout(i) = sum( vin(:) * tmp(:) ) 1750 end do 1751 end function SPC_PROJ_sp 1752 pure function SPC_PROJ_dp(space,vin) result(vout) 1753 real(dp), intent(in) :: space(:,:), vin(:) 1754 real(dp) :: vout(size(vin)), tmp(size(vin)) 1755 integer :: i 1756 do i = 1 , size(vin) 1757 tmp(:) = space(:,i) / VNORM(space(:,i)) 1758 vout(i) = sum( vin(:) * tmp(:) ) 1759 end do 1760 end function SPC_PROJ_dp 1761 1762 pure function IDX_SPC_PROJ_sp(space,vin,mag) result(idx) 1763 real(sp), intent(in) :: space(:,:), vin(:) 1764 ! If mag is false (default) it takes the largest 1765 ! POSITIVE projection. 1766 ! If mag is true it takes the largest MAGNITUDE projection 1767 logical, intent(in), optional :: mag 1768 integer :: idx 1769 real(sp) :: spc(size(vin)) 1770 integer :: i 1771 idx = 1 1772 ! Find the largest contributing space-vector 1773 ! We _only_ compare against the projected 1774 ! vector onto the normal length of the vector 1775 if ( present(mag) ) then 1776 if ( mag ) then 1777 spc(:) = abs( SPC_PROJ(space,vin) ) 1778 else 1779 spc(:) = SPC_PROJ(space,vin) 1780 end if 1781 else 1782 spc(:) = SPC_PROJ(space,vin) 1783 end if 1784 ! Find the largest contributing one 1785 ! it does not matter whether it is plus 1786 ! or minus, so long as it its magnitude 1787 ! is larger (one would _never_ have two 1788 ! vectors of opposite direction (i.e. a 2D surface) 1789 do i = 2 , size(vin) 1790 if ( spc(idx) < spc(i) ) then 1791 idx = i 1792 end if 1793 end do 1794 end function IDX_SPC_PROJ_sp 1795 1796 pure function IDX_SPC_PROJ_dp(space,vin,mag) result(idx) 1797 real(dp), intent(in) :: space(:,:), vin(:) 1798 ! If mag is false (default) it takes the largest 1799 ! POSITIVE projection. 1800 ! If mag is true it takes the largest MAGNITUDE projection 1801 logical, intent(in), optional :: mag 1802 integer :: idx 1803 real(dp) :: spc(size(vin)) 1804 integer :: i 1805 idx = 1 1806 if ( present(mag) ) then 1807 if ( mag ) then 1808 spc(:) = abs( SPC_PROJ(space,vin) ) 1809 else 1810 spc(:) = SPC_PROJ(space,vin) 1811 end if 1812 else 1813 spc(:) = SPC_PROJ(space,vin) 1814 end if 1815 do i = 2 , size(vin) 1816 if ( spc(idx) < spc(i) ) then 1817 idx = i 1818 end if 1819 end do 1820 end function IDX_SPC_PROJ_dp 1821 1822 1823 ! Scalar projection of 'vin=a' onto 'vec=b'. 1824 ! a . b / |b| 1825 ! I.e. for these vectors: 1826 ! vec = [2, 0, 0] 1827 ! vin = [3, 1, 1] 1828 ! a = [3, 0, 0] 1829 pure function VEC_PROJ_SCA_sp(vec,vin) result(a) 1830 real(sp), intent(in) :: vec(:), vin(:) 1831 real(sp) :: a 1832 a = sum(vec*vin) / VNORM(vec) 1833 end function VEC_PROJ_SCA_sp 1834 pure function VEC_PROJ_SCA_dp(vec,vin) result(a) 1835 real(dp), intent(in) :: vec(:), vin(:) 1836 real(dp) :: a 1837 a = sum(vec*vin) / VNORM(vec) 1838 end function VEC_PROJ_SCA_dp 1839 1840 ! Projection of 'vin=a' onto 'vec=b' 1841 ! a . b / ( b . b ) * b 1842 pure function VEC_PROJ_sp(vec,vin) result(vout) 1843 real(sp), intent(in) :: vec(:), vin(:) 1844 real(sp) :: vout(size(vec)) 1845 vout = sum(vec * vin) / sum(vec * vec) * vec 1846 end function VEC_PROJ_sp 1847 pure function VEC_PROJ_dp(vec,vin) result(vout) 1848 real(dp), intent(in) :: vec(:), vin(:) 1849 real(dp) :: vout(size(vec)) 1850 vout = sum(vec * vin) / sum(vec * vec) * vec 1851 end function VEC_PROJ_dp 1852 1853end module intrinsic_missing 1854 1855#ifdef INTRINSIC_MISSING_TEST 1856program test 1857 use intrinsic_missing 1858 1859 integer, parameter :: sp = selected_real_kind(5,10) 1860 integer, parameter :: dp = selected_real_kind(10,100) 1861 real(dp) :: cell(3,3), v(3) 1862 real :: t1 , t2 1863 1864 integer :: i, j 1865 integer, parameter :: N = 21355 1866 integer :: list(N) 1867 1868 1869 cell(:,:) = 0._dp 1870 cell(1,1) = 10._dp 1871 cell(1,2) = 1._dp 1872 cell(2,2) = 5._dp 1873 cell(1,3) = 1._dp 1874 cell(2,3) = 5._dp 1875 cell(3,3) = 10._dp 1876 1877 v(:) = 0._dp 1878 v(1) = 2._dp 1879 print *,1,'==',IDX_SPC_PROJ(cell,v) 1880 v(2) = 5._dp 1881 print *,2,'==',IDX_SPC_PROJ(cell,v) 1882 v(3) = 10._dp 1883 print *,3,'==',IDX_SPC_PROJ(cell,v) 1884 v(2) = v(1) 1885 v(3) = 0._dp 1886 print *,2,'==',IDX_SPC_PROJ(cell,v) 1887 v(2) = v(1) / 10._dp 1888 v(3) = 0._dp 1889 print *,1,'==',IDX_SPC_PROJ(cell,v) 1890 1891 do i = 1 , N 1892 list(i) = i 1893 end do 1894 1895 call cpu_time(t1) 1896 do j = 2 , N 1897 do i = 1 , j 1898 if ( SORTED_FIND(list(1:j),i) /= i ) print *,'Error',i 1899 end do 1900 do i = 1 , j 1901 if ( SORTED_FIND(list(1:j),i+j) /= 0 ) print *,'Error',i 1902 end do 1903 end do 1904 call cpu_time(t2) 1905 print '(a,f10.5,a)','Timing of SORTED_FIND: ',t2-t1,' secs' 1906 1907 call cpu_time(t1) 1908 do j = 2 , N 1909 do i = 1 , j 1910 if ( SORTED_REC_FIND(list(1:j),i) /= i ) print *,'Error',i 1911 end do 1912 do i = 1 , j 1913 if ( SORTED_REC_FIND(list(1:j),i+j) /= 0 ) print *,'Error',i 1914 end do 1915 end do 1916 call cpu_time(t2) 1917 print '(a,f10.5,a)','Timing of SORTED_REC_FIND: ',t2-t1,' secs' 1918 1919 call cpu_time(t1) 1920 do j = 2 , N 1921 do i = 1 , j 1922 if ( SORTED_BINARY_FIND(list(1:j),i) /= i ) print *,'Error',i 1923 end do 1924 do i = 1 , j 1925 if ( SORTED_BINARY_FIND(list(1:j),i+j) /= 0 ) print *,'Error',i 1926 end do 1927 end do 1928 call cpu_time(t2) 1929 print '(a,f10.5,a)','Timing of SORTED_BINARY_FIND: ',t2-t1,' secs' 1930 1931end program test 1932#endif 1933