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