1! $Id: restart.f90,v 1.13 2007/06/25 20:16:31 kewu Exp $ 2!!! TRLAN Low level utility routines 3! This file contains a number of routines related to decisions of how 4! many Ritz pairs to save when restarting the Thick-Restart Lanczos 5! method. The subroutine trl_shuffle_eig is the main access point. 6!!! 7! The subroutine trl_shuffle_eig accomplishes two tasks: 8! (1) sort the Ritz values so that those to be saved after 9! restart are ordered in the front of the array, 10! (2) determine the number of Ritz values to be saved. 11! On return from this subroutine, the Ritz values are in ascending 12! order so that DSTEIN can be used to compute the eigenvectors 13!!! 14Subroutine trl_shuffle_eig(nd, mnd, lambda, res, info, kept) 15 Use trl_info 16 Implicit None 17 Type(TRL_INFO_T), Intent(in) :: info 18 Integer, Intent(in) :: nd, mnd 19 Integer :: kept 20 Real(8) :: lambda(nd), res(nd) 21 ! 22 Integer :: i, ncl, ncr, kl, kr, tind, minsep 23 Real(8) :: bnd 24 External dsort2 25 Interface 26 Subroutine trl_restart_fixed(nd, mnd, tind, lambda, res, info, kl, kr) 27 Use trl_info 28 Implicit None 29 Type(TRL_INFO_T), Intent(in) :: info 30 Integer, Intent(in) :: nd, mnd, tind 31 Integer, Intent(inout) :: kl, kr 32 Real(8), Intent(in) :: lambda(nd), res(nd) 33 End Subroutine trl_restart_fixed 34 Subroutine trl_restart_scan(nd, res, info, kept, kl, kr) 35 Use trl_info 36 Implicit None 37 Type(TRL_INFO_T), Intent(in) :: info 38 Integer, Intent(in) :: nd, kept 39 Integer, Intent(inout) :: kl, kr 40 Real(8), Intent(in) :: res(nd) 41 End Subroutine trl_restart_scan 42 Subroutine trl_restart_small_res(nd, mnd, tind, lambda, res, info, kl, kr) 43 Use trl_info 44 Implicit None 45 Type(TRL_INFO_T), Intent(in) :: info 46 Integer, Intent(in) :: nd, mnd, tind 47 Integer, Intent(inout) :: kl, kr 48 Real(8), Intent(in) :: lambda(nd), res(nd) 49 End Subroutine trl_restart_small_res 50 Subroutine trl_restart_max_gap_ratio(nd, tind, kept, lambda, res,& 51 & info, kl, kr) 52 Use trl_info 53 Implicit None 54 Type(TRL_INFO_T), Intent(in) :: info 55 Integer, Intent(in) :: kept, nd, tind 56 Integer, Intent(inout) :: kl, kr 57 Real(8), Intent(in) :: lambda(nd), res(nd) 58 End Subroutine trl_restart_max_gap_ratio 59 Subroutine trl_restart_max_progress(nd, tind, kept, lambda, res,& 60 & info, kl, kr) 61 Use trl_info 62 Implicit None 63 Type(TRL_INFO_T), Intent(in) :: info 64 Integer, Intent(in) :: kept, nd, tind 65 Integer, Intent(inout) :: kl, kr 66 Real(8), Intent(in) :: lambda(nd), res(nd) 67 End Subroutine trl_restart_max_progress 68 Subroutine trl_restart_max_reduction(nd, tind, kept, lambda, res,& 69 & info, kl, kr) 70 Use trl_info 71 Implicit None 72 Type(TRL_INFO_T), Intent(in) :: info 73 Integer, Intent(in) :: kept, nd, tind 74 Integer, Intent(inout) :: kl, kr 75 Real(8), Intent(in) :: lambda(nd), res(nd) 76 End Subroutine trl_restart_max_reduction 77 End Interface 78 ! 79 ! very small basis -- save the half with the smallest residual norms 80 If (nd .Le. 5) Then 81 Call dsort2(nd, res, lambda) 82 If (nd.Gt.3) Then 83 kept = 2 84 Else If (nd.Gt.0) Then 85 kept = 1 86 Else 87 kept = 0 88 End If 89 If (kept.Gt.1) Call dsort2(kept, lambda, res) 90 Return 91 End If 92 ! 93 ! preparation for normal case, first determine what are converged 94 ! ncl are converged from the left and ncr are converged from the 95 ! right 96 bnd = Min(info%tol, Epsilon(info%tol))*info%anrm 97 ncr = 1 98 ncl = nd 99 i = nd 100 Do While (i .Gt. 0) ! determine how many has converged from the right 101 If (res(i) .Le. bnd) Then 102 i = i - 1 103 Else 104 ncr = i + 1 105 i = 0 106 End If 107 End Do 108 i = 1 109 Do While (i .Le. nd) ! determine how many has converged from the left 110 If (res(i) .Le. bnd) Then 111 i = i + 1 112 Else 113 ncl = i - 1 114 i = nd + 1 115 End If 116 End Do 117 kl = Max(1, ncl) 118 kr = Min(nd, ncr) 119 If (ncr .Gt. ncl) Then 120 ! find the one that is closest to info%trgt 121 tind = (kl+kr)/2 122 Do While (lambda(tind).Ne.info%trgt .And. kr.Gt.kl) 123 If (lambda(tind) .Lt. info%trgt) Then 124 kl = tind + 1 125 tind = (kl + kr) / 2 126 Else If (lambda(tind) .Gt. info%trgt) Then 127 kr = tind - 1 128 tind = (kl + kr) / 2 129 Else 130 kl = tind 131 kr = tind 132 End If 133 End Do 134 ! assign kl to the largest index of lambda that is smaller than 135 ! info%trgt 136 If (lambda(tind) .Eq. info%trgt) Then 137 kl = tind - 1 13810 If (kl .Gt. 0) Then 139 If (lambda(kl) .Eq. info%trgt) Then 140 kl = kl - 1 141 Goto 10 142 End If 143 End If 144 ! assign kr to the smallest index of lambda that is greater than 145 ! info%trgt 146 kr = tind + 1 14720 If (kr .Le. nd) Then 148 If (lambda(kr) .Eq. info%trgt) Then 149 kr = kr + 1 150 Goto 20 151 End If 152 End If 153 Else 154 kl = tind - 1 155 kr = tind + 1 156 End If 157 ! initial assignment of kl and kr 158 If (info%lohi .Gt. 0) Then 159 kr = kl 160 kl = Min(ncl, Max(0, nd-info%ned)) 161 Else If (info%lohi .Lt. 0) Then 162 kl = kr 163 kr = Max(ncr, Min(nd-info%nec,info%ned+1)) 164 Else If (ncr-tind .Gt. tind-ncl) Then 165 kl = kr 166 kr = Max(ncr, Min(nd-info%nec,info%ned+1)) 167 Else 168 kr = kl 169 kl = Min(ncl, Max(0, nd-info%ned)) 170 End If 171 Else 172 ! all have converged, keep all -- should not happen 173 kept = nd 174 Return 175 End If 176 ! 177 ! normal cases, call subroutines to complete the tasks 178 ! [1 .. kl] and [kr .. nd] are saved for later 179 ! the initial values of kl and kr are simply ncl and ncr 180 ! they are further analyzed according to the restarting strategy 181 ! requested 182 Select Case (info%restart) 183 Case (1) 184 ! fixed number beyond the currently converged ones 185 Call trl_restart_fixed(nd, mnd, tind, lambda, res, info, kl, kr) 186 Case (2) 187 ! add the ones with smallest residual nroms 188 Call trl_restart_small_res(nd, mnd, tind, lambda, res, info, kl, kr) 189 Case (3) 190 If (info%nloop .Gt. 0) Then 191 ! maximize the gap ratio 192 Call trl_restart_max_gap_ratio(nd, tind, kept, lambda, res, info,& 193 & kl, kr) 194 Else 195 ! this is the first restart -- use trl_restart_fixed instead 196 Call trl_restart_fixed(nd, mnd, tind, lambda, res, info, kl, kr) 197 End If 198 Case (4) 199 If (info%nloop .Gt. 0) Then 200 ! maximize [gap-ratio * (m-k)] 201 Call trl_restart_max_progress(nd, tind, kept, lambda, res,& 202 & info, kl, kr) 203 Else 204 ! this is the first restart -- use trl_restart_fixed instead 205 Call trl_restart_fixed(nd, mnd, tind, lambda, res, info, kl, kr) 206 End If 207 Case (5) 208 If (info%nloop .Gt. 0) Then 209 ! maximize [sqrt(gap tatio) * (m-k)] 210 Call trl_restart_max_reduction(nd, tind, kept, lambda, res,& 211 & info, kl, kr) 212 Else 213 ! this is the first restart -- use trl_restart_fixed instead 214 Call trl_restart_fixed(nd, mnd, tind, lambda, res, info, kl, kr) 215 End If 216 Case (6) 217 ! progressively vary the thickness 218 Call trl_restart_scan(nd, res, info, kept, kl, kr) 219 Case default 220 If (info%restart .Le. -info%ned) Then 221 If (info%lohi .Ge. 0) Then 222 kl = 0 223 kr = Max(2, nd+info%restart)+1 224 Else If (info%lohi .Lt. 0) Then 225 kl = Min(-info%restart, nd-3) 226 kr = nd+1 227 Else 228 kl = Min(nd-3, -info%restart)/2 229 kr = nd - kl 230 End If 231 Else 232 Call trl_restart_fixed(nd, mnd, tind, lambda, res, info, kl, kr) 233 End If 234 End Select 235 ! 236 ! make sure kr > kl+minsep 237 minsep = Max(3, nd/6, nd-6*info%ned) 238 If (kr.Le.kl+minsep .Or. kl+nd-kr+minsep.Gt.mnd) Then 239 If (ncl.Lt.kl .And. kl.Lt.kr .And. kr.Lt.ncr) Then 240 kl = kl - 1 241 kr = kr + 1 242 Else If (info%lohi .Gt. 0) Then 243 kr = Max(minsep, Min(nd/3, ncr-1)) 244 kl = 0 245 Else If (info%lohi .Lt. 0) Then 246 kl = Min(nd-minsep, Max((nd+nd)/3, ncl+1)) 247 kr = nd+1 248 Else 249 kl = (nd-minsep)/2-1 250 kr = (nd-minsep+1)/2+1 251 End If 252 End If 253 ! copy the (kr:nd) elements to (kl+1:kl+nd-kr+2) 254 ! kr is temporarily decreased by 1 to make indices easier to compute 255 kr = kr - 1 256 Do i = 1, nd-kr 257 lambda(kl+i) = lambda(kr+i) 258 res(kl+i) = res(kr+i) 259 End Do 260 kept = kl + Max(0, nd-kr) 261 Return 262End Subroutine trl_shuffle_eig 263!!! 264! save fixed number of extra Ritz pairs 265! January 99 -- added modification to ensure a minimal gap ratio is 266! maintained 267!!! 268Subroutine trl_restart_fixed(nd, mnd, tind, lambda, res, info, kl, kr) 269 Use trl_info 270 Implicit None 271 Type(TRL_INFO_T), Intent(in) :: info 272 Integer, Intent(in) :: nd, mnd, tind 273 Integer, Intent(inout) :: kl, kr 274 Real(8), Intent(in) :: lambda(nd), res(nd) 275 Interface 276 Function trl_min_gap_ratio(info, nd, tind, res) Result(gamma) 277 Use trl_info 278 Implicit None 279 Type(TRL_INFO_T), Intent(in) :: info 280 Integer, Intent(in) :: nd, tind 281 Real(8), Intent(in) :: res(nd) 282 Real(8) :: gamma 283 End Function trl_min_gap_ratio 284 End Interface 285 ! 286 ! local variables 287 Real(8), Parameter :: ten=1.0D1 288 Integer :: extra, i, kl0, kr0, minsep 289 Real(8) :: gmin 290 ! the number of extra Ritz pairs to be saved 291 kl0 = kl 292 kr0 = kr 293 extra = Nint((mnd-info%nec)*(0.4D0+0.1D0*info%ned/Dble(mnd))) 294 If (extra.Gt.info%ned+info%ned .And. extra.Gt.5) Then 295 gmin = Dble(mnd)/Dble(info%ned) 296 extra = Nint((extra + (Log(gmin)*info%ned) * gmin) / (1.0D0 +& 297 & gmin)) 298 End If 299 minsep = Max(3, nd/5, nd-4*info%ned) 300 gmin = trl_min_gap_ratio(info, nd, tind, res) 301 If (info%lohi .Gt. 0) Then 302 kr = Min(tind-1, kr) - extra 303 kl = 0 304 Else If (info%lohi .Lt. 0) Then 305 kl = Max(tind+1, kl) + extra 306 kr = nd+1 307 Else 308 extra = extra + 1 309 kl = kl + extra/2 310 kr = kr - extra/2 311 i = 1 312 Do While (kl.Gt.kl0 .And. kr.Lt.kr0 .And. i.Eq.1) 313 If (ten*res(kl).Lt.res(kr)) Then 314 If (res(kl+1).Lt.res(kr+1)) Then 315 kl = kl + 1 316 kr = kr + 1 317 Else 318 i = 0 319 End If 320 Else If (ten*res(kr).Lt.res(kl)) Then 321 If (res(kr-1).Lt.res(kl-1)) Then 322 kr = kr - 1 323 kl = kl - 1 324 Else 325 i = 0 326 End If 327 Else 328 i = 0 329 End If 330 End Do 331 End If 332 ! adjust kl and kr until the minimal gap ratio is satisfied 333 Do While (kl+minsep.Lt.kr .And. gap_ratio(Max(1,kl),Min(kr,nd)).Lt.gmin) 334 If (info%lohi .Lt. 0) Then 335 kl = kl + 1 336 Else If (info%lohi .Gt. 0) Then 337 kr = kr - 1 338 Else If (res(kl) .Lt. res(kr)) Then 339 kl = kl + 1 340 Else 341 kr = kr + 1 342 End If 343 End Do 344 ! make sure nearly degenerate Ritz pairs are included 345 ! lambda(kl)+r(kl) > lambda(j) and 346 ! lambda(kl)-r(kl) > lambda(j)-r(j) 347 If (info%lohi .Gt. 0) Then 348 i = kr-1 349 Do While (i.Gt.kl+minsep .And. lambda(kr)-res(kr).Lt.lambda(i) .And.& 350 & lambda(kr)+res(kr).Lt.lambda(i)+res(i)) 351 i = i - 1 352 End Do 353 kr = i+1 354 Else 355 kl0 = kl 356 i = kl + 1 357 Do While (i.Lt.kr-minsep .And. lambda(kl)+res(kl).Gt.lambda(i) .And.& 358 & lambda(kl)-res(kl).Gt.lambda(i)-res(i)) 359 i = i + 1 360 End Do 361 kl = i-1 362 End If 363Contains 364 Function gap_ratio(i,j) Result(gamma) 365 Integer, Intent(in) :: i, j 366 Real(8) :: gamma 367 gamma = (lambda(i) - lambda(tind)) / (lambda(j) - lambda(tind)) 368 Return 369 End Function gap_ratio 370End Subroutine trl_restart_fixed 371!!! 372! progressively vary the thickness to scan all possible values 373! 374! This subroutine determines the number Ritz vectors to be saved by 375! adding some quantity to the current thickness. If the thickness is 376! larger than nd-2, it is reset to something smaller. 377!!! 378Subroutine trl_restart_scan(nd, res, info, kept, kl, kr) 379 Use trl_info 380 Implicit None 381 Type(TRL_INFO_T), Intent(in) :: info 382 Integer, Intent(in) :: nd, kept 383 Integer, Intent(inout) :: kl, kr 384 Real(8), Intent(in) :: res(nd) 385 ! 386 ! local variables 387 Real(8), Parameter :: ten=1.0D1 388 Integer :: extra, i, kl0, kr0 389 ! three cases have to be dealt with separately 390 If (info%lohi .Lt. 0) Then 391 kr = nd + 1 392 kl = kept + Min(Max(info%nec,1), (nd-kept)/2) 393 If (kl .Le. 1) Then 394 If (nd.Gt.6) Then 395 kl = nd/2 396 Else If (nd.Gt.2) Then 397 kl = 2 398 End If 399 Else If (kl+3 .Ge. nd) Then 400 kl = info%nec + Min(info%ned, 10, (nd-info%ned)/2) 401 End If 402 Else If (info%lohi .Gt. 0) Then 403 kl = 0 404 kr = kept + Min(Max(info%nec,1), (nd-kept)/2) 405 If (kr .Le. 1) Then 406 If (nd .Gt. 6) Then 407 kr = nd / 2 408 Else If (nd .Gt. 2) Then 409 kr = 2 410 End If 411 Else If (kr+3 .Gt. nd) Then 412 kr = info%nec + Min(info%ned, 10, (nd-info%ned)/2) 413 End If 414 kr = nd - kr + 1 415 Else 416 kl0 = kl 417 kr0 = kr 418 extra = kept + Min(info%nec, (nd-kept)/2) + 1 419 If (extra .Le. 1) Then 420 If (nd .Gt. 6) Then 421 extra = nd / 2 422 Else If (nd .Gt. 2) Then 423 extra = 2 424 End If 425 Else If (extra+3 .Gt. nd) Then 426 extra = info%nec + Min(info%ned, 10, (nd-info%ned)/2) 427 End If 428 kl = Max(kl, extra/2) 429 kr = Min(kr, nd+1-extra/2) 430 i = 1 431 Do While (kl.Gt.kl0 .And. kr.Lt.kr0 .And. i.Eq.1) 432 If (ten*res(kl).Lt.res(kr)) Then 433 If (res(kl+1).Lt.res(kr+1)) Then 434 kl = kl + 1 435 kr = kr + 1 436 Else 437 i = 0 438 End If 439 Else If (ten*res(kr).Lt.res(kl)) Then 440 If (res(kr-1).Lt.res(kl-1)) Then 441 kr = kr - 1 442 kl = kl - 1 443 Else 444 i = 0 445 End If 446 Else 447 i = 0 448 End If 449 End Do 450 End If 451End Subroutine trl_restart_scan 452!!! 453! save those that are close to converge (as close as the target) 454!!! 455Subroutine trl_restart_small_res(nd, mnd, tind, lambda, res, info, kl, kr) 456 Use trl_info 457 Implicit None 458 Type(TRL_INFO_T), Intent(in) :: info 459 Integer, Intent(in) :: nd, mnd, tind 460 Integer, Intent(inout) :: kl, kr 461 Real(8), Intent(in) :: lambda(nd), res(nd) 462 Interface 463 Function trl_min_gap_ratio(info, nd, tind, res) Result(gamma) 464 Use trl_info 465 Implicit None 466 Type(TRL_INFO_T), Intent(in) :: info 467 Integer, Intent(in) :: nd, tind 468 Real(8), Intent(in) :: res(nd) 469 Real(8) :: gamma 470 End Function trl_min_gap_ratio 471 End Interface 472 ! 473 ! local variables 474 Integer :: extra, i, j, ii, kl0, kr0, minsep 475 Real(8) :: acpt, resmax, gmin 476 ! the number of extra Ritz pairs to be saved 477 minsep = Max(3, nd/5, nd-4*info%ned) 478 extra = Nint((mnd-info%nec)*(0.4D0+0.1D0*info%ned/Dble(mnd))) 479 If (extra.Gt.info%ned+info%ned .And. extra.Gt.5) Then 480 gmin = Dble(mnd)/Dble(info%ned) 481 extra = Nint((extra + (Log(gmin)*info%ned) * gmin) / (1.0D0 +& 482 & gmin)) 483 End If 484 kl0 = kl 485 kr0 = kr 486 resmax = Maxval(res) 487 acpt = resmax / res(tind) 488 ! 489 ! determine the number of Ritz pairs that has to be saved 490 If (info%lohi .Gt. 0) Then 491 If (acpt.Lt.0.999D0 .And. acpt.Ge.0.0D0) Then 492 ii = tind - 1 493 acpt = Max(Sqrt(acpt)*res(tind), res(ii)+res(ii)) 494 acpt = Min(acpt, resmax) 495 kr = ii - 1 496 Do While (res(kr).Lt.acpt .And. kr.Gt.kl+3) 497 kr = kr - 1 498 End Do 499 Else 500 kr = kr0 - extra 501 End If 502 kr = Max(3, kr) 503 kl = Min(kl, kr-2) 504 Else If (info%lohi .Lt. 0) Then 505 If (acpt.Lt.0.999D0 .And. acpt.Ge.0.0D0) Then 506 ii = tind + 1 507 acpt = Max(Sqrt(acpt)*res(tind), res(ii)+res(ii)) 508 acpt = Min(acpt, resmax) 509 kl = ii + 1 510 Do While (res(kl).Lt.acpt .And. kl.Lt.kr-3) 511 kl = kl + 1 512 End Do 513 Else 514 kl = kl0 + extra 515 End If 516 kl = Min(nd-3, kl) 517 kr = Max(kr, kl+2) 518 Else 519 ! save whoever has smaller residual norms 520 i = kl + 1 521 j = kr - 1 522 Do ii = 1, extra 523 If (res(i) .Lt. res(j)) Then 524 kl = i 525 i = i + 1 526 Else If (res(i) .Gt. res(j)) Then 527 kr = j 528 j = j - 1 529 Else If (i .Le. nd-j) Then 530 kl = i 531 i = i + 1 532 Else 533 kr = j 534 j = j - 1 535 End If 536 End Do 537 End If 538 ! adjust kl and kr until the minimal gap ratio is satisfied 539 kl0 = kl 540 kr0 = kr 541 gmin = trl_min_gap_ratio(info, nd, tind, res) 542 Do While (kl+minsep.Lt.kr .And. gap_ratio(Max(1,kl),Min(nd,kr)).Lt.gmin) 543 If (info%lohi .Lt. 0) Then 544 kl = kl + 1 545 Else If (info%lohi .Gt. 0) Then 546 kr = kr - 1 547 Else If (res(kl) .Lt. res(kr)) Then 548 kl = kl + 1 549 Else 550 kr = kr + 1 551 End If 552 End Do 553 ! make sure nearly degenerate Ritz pairs are included 554 ! lambda(kl)+r(kl) > lambda(j) and 555 ! lambda(kl)-r(kl) > lambda(j)-r(j) 556 If (info%lohi .Gt. 0) Then 557 i = kr0-1 558 Do While (i.Gt.kl+minsep .And. lambda(kr)-res(kr).Lt.lambda(i) .And.& 559 & lambda(kr)+res(kr).Lt.lambda(i)+res(i)) 560 i = i - 1 561 End Do 562 kr = Min(kr, i+1) 563 Else 564 i = kl0 + 1 565 Do While (i.Lt.kr-minsep .And. lambda(kl)+res(kl).Gt.lambda(i) .And.& 566 & lambda(kl)-res(kl).Gt.lambda(i)-res(i)) 567 i = i + 1 568 End Do 569 kl = Max(kl, i-1) 570 End If 571Contains 572 Function gap_ratio(i,j) Result(gamma) 573 Integer, Intent(in) :: i, j 574 Real(8) :: gamma 575 gamma = (lambda(i) - lambda(tind)) / (lambda(j) - lambda(tind)) 576 Return 577 End Function gap_ratio 578End Subroutine trl_restart_small_res 579!!! 580! search throught all pairs of (kl, kr) for the one with the maximum 581! gap ratio for the next Ritz pair(target) 582! 583! This is an optimized version of the original version. It only search 584! through nd number once. (Only single loop!) 585!!! 586Subroutine trl_restart_max_gap_ratio(nd, tind, kept, lambda, res, info, kl, kr) 587 Use trl_info 588 Implicit None 589 Type(TRL_INFO_T), Intent(in) :: info 590 Integer, Intent(in) :: kept, nd, tind 591 Integer, Intent(inout) :: kl, kr 592 Real(8), Intent(in) :: lambda(nd), res(nd) 593 Interface 594 Subroutine trl_restart_search_range(nd, lambda, res, info, ncl,& 595 & ncr, lohi, tind, klm, krm) 596 Use trl_info 597 Implicit None 598 Type(TRL_INFO_T), Intent(in) :: info 599 Integer, Intent(in) :: nd, ncl, ncr, tind 600 Integer, Intent(out) :: klm, krm, lohi 601 Real(8), Intent(in) :: lambda(nd), res(nd) 602 End Subroutine trl_restart_search_range 603 Function trl_min_gap_ratio(info, nd, tind, res) Result(gamma) 604 Use trl_info 605 Implicit None 606 Type(TRL_INFO_T), Intent(in) :: info 607 Integer, Intent(in) :: nd, tind 608 Real(8), Intent(in) :: res(nd) 609 Real(8) :: gamma 610 End Function trl_min_gap_ratio 611 End Interface 612 ! 613 ! local variables 614 Integer :: i, j, lohi, klm, krm, igap 615 Real(8) :: bnd, tmp, gratio 616 ! statement function for computing gap ratio 617 gratio(i,j) = (lambda(i)-info%trgt)/(lambda(j)-info%trgt) 618 ! determine the search range 619 Call trl_restart_search_range(nd, lambda, res, info, kl, kr, lohi,& 620 & tind, klm, krm) 621 kl = klm 622 kr = krm 623 igap = Max(Min(nd-info%ned, Nint((krm-klm)*0.4D0)), 2) 624 If (igap.Gt.2 .And. info%matvec.Gt.info%maxlan) Then 625 If (info%clk_op+info%tick_o.Gt.1.0D1*(info%clk_orth& 626 &+info%tick_h+info%clk_res+info%tick_r)) Then 627 igap = Max(2, nd-kept-1) 628 Else 629 bnd = trl_min_gap_ratio(info, nd, tind, res) 630 If (info%crate .Lt. bnd) igap = Max(2, nd-kept-1) 631 End If 632 End If 633 If (kl+igap .Gt. kr) Return 634 ! two cases depending on lohi 635 If (lohi .Gt. 0) Then 636 ! target is at the high end of spectrum 637 bnd = gratio(kr, kl) 638 Do i = klm, krm-igap 639 j = i + igap 640 tmp = gratio(j,i) 641 If (tmp .Gt. bnd) Then 642 kl = i 643 kr = j 644 bnd = tmp 645 End If 646 End Do 647 Else 648 bnd = gratio(kl, kr) 649 Do i = klm, krm-igap 650 j = i + igap 651 tmp = gratio(i,j) 652 If (tmp .Gt. bnd) Then 653 kl = i 654 kr = j 655 bnd = tmp 656 End If 657 End Do 658 End If 659End Subroutine trl_restart_max_gap_ratio 660!!! 661! search for a pair (kl, kr) such that the reduction in residual norm 662! of the target (info%trgt) will be the largest before next restart 663! The merit function is [gap-ratio * (m-k)] 664!!! 665Subroutine trl_restart_max_progress(nd, tind, kept, lambda, res, info, kl, kr) 666 Use trl_info 667 Implicit None 668 Type(TRL_INFO_T), Intent(in) :: info 669 Integer, Intent(in) :: kept, nd, tind 670 Integer, Intent(inout) :: kl, kr 671 Real(8), Intent(in) :: lambda(nd), res(nd) 672 Interface 673 Subroutine trl_restart_search_range(nd, lambda, res, info, ncl,& 674 & ncr, lohi, tind, klm, krm) 675 Use trl_info 676 Implicit None 677 Type(TRL_INFO_T), Intent(in) :: info 678 Integer, Intent(in) :: nd, ncl, ncr, tind 679 Integer, Intent(out) :: klm, krm, lohi 680 Real(8), Intent(in) :: lambda(nd), res(nd) 681 End Subroutine trl_restart_search_range 682 Function trl_min_gap_ratio(info, nd, tind, res) Result(gamma) 683 Use trl_info 684 Implicit None 685 Type(TRL_INFO_T), Intent(in) :: info 686 Integer, Intent(in) :: nd, tind 687 Real(8), Intent(in) :: res(nd) 688 Real(8) :: gamma 689 End Function trl_min_gap_ratio 690 End Interface 691 ! 692 ! local variables 693 Integer :: i, j, lohi, klm, krm, igap 694 Real(8) :: tmp, ss, merit 695 ! the merit function for determining how to restart 696 merit(i,j) = (lambda(i)-info%trgt) * Abs(j-i) / (lambda(j)-info%trgt) 697 ! determine the search range 698 Call trl_restart_search_range(nd, lambda, res, info, kl, kr, lohi,& 699 & tind, klm, krm) 700 ! perform the brute-force search 701 kl = klm 702 kr = krm 703 igap = Max(Min(nd-info%ned, Nint((krm-klm)*0.4D0)), 2) 704 If (igap.Gt.2 .And. igap+kept.Gt.nd .And. info%crate.Gt.0.0D0) Then 705 ss = trl_min_gap_ratio(info, nd, tind, res) 706 If (ss .Gt. info%crate) igap = Max(2, nd-kept-1) 707 End If 708 If (lohi .Gt. 0) Then 709 ss = merit(kr, kl) 710 Do i = klm, krm-igap 711 Do j = i+igap, krm 712 tmp = merit(j, i) 713 If (tmp.Gt.ss) Then 714 ss = tmp 715 kl = i 716 kr = j 717 End If 718 End Do 719 End Do 720 Else 721 ss = merit(kl, kr) 722 Do i = klm, krm-igap 723 Do j = i+igap, krm 724 tmp = merit(i,j) 725 If (tmp.Gt.ss) Then 726 ss = tmp 727 kl = i 728 kr = j 729 End If 730 End Do 731 End Do 732 End If 733End Subroutine trl_restart_max_progress 734!!! 735! search for a pair (kl, kr) such that the reduction in residual norm 736! of the target (info%trgt) will be the largest before next restart 737! the merit function is [sqrt(gap ratio) * (m-k)] 738!!! 739Subroutine trl_restart_max_reduction(nd, tind, kept, lambda, res, info, kl, kr) 740 Use trl_info 741 Implicit None 742 Type(TRL_INFO_T), Intent(in) :: info 743 Integer, Intent(in) :: kept, nd, tind 744 Integer, Intent(inout) :: kl, kr 745 Real(8), Intent(in) :: lambda(nd), res(nd) 746 Interface 747 Subroutine trl_restart_search_range(nd, lambda, res, info, ncl,& 748 & ncr, lohi, tind, klm, krm) 749 Use trl_info 750 Implicit None 751 Type(TRL_INFO_T), Intent(in) :: info 752 Integer, Intent(in) :: nd, ncl, ncr, tind 753 Integer, Intent(out) :: klm, krm, lohi 754 Real(8), Intent(in) :: lambda(nd), res(nd) 755 End Subroutine trl_restart_search_range 756 Function trl_min_gap_ratio(info, nd, tind, res) Result(gamma) 757 Use trl_info 758 Implicit None 759 Type(TRL_INFO_T), Intent(in) :: info 760 Integer, Intent(in) :: nd, tind 761 Real(8), Intent(in) :: res(nd) 762 Real(8) :: gamma 763 End Function trl_min_gap_ratio 764 End Interface 765 ! 766 ! local variables 767 Integer :: i, j, lohi, klm, krm, igap 768 Real(8) :: tmp, ss, merit 769 ! merit measure the factor of residual norm reduction 770 merit(i,j) = Sqrt((lambda(i)-info%trgt)/(lambda(j)-info%trgt)) * Abs(j-i) 771 ! determine the search range 772 Call trl_restart_search_range(nd, lambda, res, info, kl, kr, lohi,& 773 & tind, klm, krm) 774 ! perform the brute-force search 775 kl = klm 776 kr = krm 777 igap = Max(Min(nd-info%ned, Nint((krm-klm)*0.4D0)), 2) 778 If (igap.Gt.2 .And. igap+kept.Gt.nd .And. info%crate.Gt.0.0D0) Then 779 ss = trl_min_gap_ratio(info, nd, tind, res) 780 If (ss .Gt. info%crate) igap = Max(2, nd-kept-1) 781 End If 782 If (lohi .Gt. 0) Then 783 ss = merit(kr, kl) 784 Do i = klm, krm-igap 785 Do j = i+igap, krm 786 tmp = merit(j, i) 787 If (tmp.Gt.ss) Then 788 ss = tmp 789 kl = i 790 kr = j 791 End If 792 End Do 793 End Do 794 Else 795 ss = merit(kl, kr) 796 Do i = klm, krm-igap 797 Do j = i+igap, krm 798 tmp = merit(i,j) 799 If (tmp.Gt.ss) Then 800 ss = tmp 801 kl = i 802 kr = j 803 End If 804 End Do 805 End Do 806 End If 807End Subroutine trl_restart_max_reduction 808!!! 809! determine the search range -- used by the schemes that performs brute 810! -force search 811!!! 812Subroutine trl_restart_search_range(nd, lambda, res, info, ncl, ncr,& 813 & lohi, tind, klm, krm) 814 Use trl_info 815 Implicit None 816 Type(TRL_INFO_T), Intent(in) :: info 817 Integer, Intent(in) :: nd, ncl, ncr, tind 818 Integer, Intent(out) :: klm, krm, lohi 819 Real(8), Intent(in) :: lambda(nd), res(nd) 820 ! 821 ! local variables 822 Integer :: j 823 Real(8) :: bnd 824 klm = Max(ncl,1) 825 krm = Min(ncr,nd) 826 bnd = info%tol * info%anrm 827 lohi = info%lohi 828 ! make sure there is some distance between the boundary and the 829 ! target Ritz value 830 If (info%lohi .Gt. 0) Then 831 ! save high end 832 krm = Min(Max(info%maxlan-info%ned, (info%maxlan+info%nec)/2),& 833 & krm, tind-1) 834 Do While (krm+krm.Ge.ncl+ncr .And. res(krm).Le.bnd) 835 krm = krm - 1 836 End Do 837 Else If (info%lohi .Lt. 0) Then 838 ! save lower end 839 klm = Max(Min(info%ned, (info%maxlan+info%nec)/2), tind+1, klm) 840 Do While (klm+klm.Le.ncl+ncr .And. res(klm).Le.bnd) 841 klm = klm + 1 842 End Do 843 Else 844 ! save both ends 845 If (tind-klm .Lt. krm-tind) Then 846 lohi = -1 847 klm = tind + 1 848 Else 849 lohi = 1 850 krm = tind - 1 851 End If 852 j = info%locked + klm + nd - krm + 1 853 If (j .Gt. 0) Then 854 j = j / 2 855 klm = klm - j 856 krm = krm + j 857 End If 858 End If 859 If (klm .Lt. 1) klm = 1 860 If (krm .Gt. nd) krm = nd 861End Subroutine trl_restart_search_range 862!!! 863! determine the minimal gap ratio needed to compute all wanted 864! eigenvalues 865!!! 866Function trl_min_gap_ratio(info, nd, tind, res) Result(gamma) 867 Use trl_info 868 Implicit None 869 Type(TRL_INFO_T), Intent(in) :: info 870 Integer, Intent(in) :: nd, tind 871 Real(8), Intent(in) :: res(nd) 872 Real(8) :: gamma 873 ! 874 gamma = info%maxmv * (info%nec + 1.0D0) / info%ned - info%matvec 875 If (gamma .Lt. info%klan) Then 876 gamma = Max((info%maxmv-info%matvec)/Dble(info%ned-info%nec), 2.0D0) 877 End If 878 gamma = Min(Log(res(tind)/(info%tol*info%anrm))/gamma, 0.5D0) 879 Return 880End Function trl_min_gap_ratio 881