restart.f90

来自「用于计算矩阵的特征值,及矩阵的其他运算.可以用与稀疏矩阵」· F90 代码 · 共 879 行 · 第 1/2 页

              kr = kr - 1
              kl = kl - 1
           Else
              i = 0
           End If
        Else
           i = 0
        End If
     End Do
  End If
End Subroutine trl_restart_scan
!!!
! save those that are close to converge (as close as the target)
!!!
Subroutine trl_restart_small_res(nd, mnd, tind, lambda, res, info, kl, kr)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T), Intent(in) :: info
  Integer, Intent(in) :: nd, mnd, tind
  Integer, Intent(inout) :: kl, kr
  Real(8), Intent(in) :: lambda(nd), res(nd)
  Interface
     Function trl_min_gap_ratio(info, nd, tind, res) Result(gamma)
       Use trl_info
       Implicit None
       Type(TRL_INFO_T), Intent(in) :: info
       Integer, Intent(in) :: nd, tind
       Real(8), Intent(in) :: res(nd)
       Real(8) :: gamma
     End Function trl_min_gap_ratio
  End Interface
  !
  ! local variables
  Integer :: extra, i, j, ii, kl0, kr0, minsep
  Real(8) :: acpt, resmax, gmin
  ! the number of extra Ritz pairs to be saved
  minsep = Max(3, nd/5, nd-4*info%ned)
  extra = Nint((mnd-info%nec)*(0.4D0+0.1D0*info%ned/Dble(mnd)))
  If (extra.Gt.info%ned+info%ned .And. extra.Gt.5) Then
     gmin = Dble(mnd)/Dble(info%ned)
     extra = Nint((extra + (Log(gmin)*info%ned) * gmin) / (1.0D0 +&
          & gmin))
  End If
  kl0 = kl
  kr0 = kr
  resmax = Maxval(res)
  acpt = resmax / res(tind)
  !
  ! determine the number of Ritz pairs that has to be saved
  If (info%lohi .Gt. 0) Then
     If (acpt.Lt.0.999D0 .And. acpt.Ge.0.0D0) Then
        ii = tind - 1
        acpt = Max(Sqrt(acpt)*res(tind), res(ii)+res(ii))
        acpt = Min(acpt, resmax)
        kr = ii - 1
        Do While (res(kr).Lt.acpt .And. kr.Gt.kl+3)
           kr = kr - 1
        End Do
     Else
        kr = kr0 - extra
     End If
     kr = Max(3, kr)
     kl = Min(kl, kr-2)
  Else If (info%lohi .Lt. 0) Then
     If (acpt.Lt.0.999D0 .And. acpt.Ge.0.0D0) Then
        ii = tind + 1
        acpt = Max(Sqrt(acpt)*res(tind), res(ii)+res(ii))
        acpt = Min(acpt, resmax)
        kl = ii + 1
        Do While (res(kl).Lt.acpt .And. kl.Lt.kr-3)
           kl = kl + 1
        End Do
     Else
        kl = kl0 + extra
     End If
     kl = Min(nd-3, kl)
     kr = Max(kr, kl+2)
  Else
     ! save whoever has smaller residual norms
     i = kl + 1
     j = kr - 1
     Do ii = 1, extra
        If (res(i) .Lt. res(j)) Then
           kl = i
           i = i + 1
        Else If (res(i) .Gt. res(j)) Then
           kr = j
           j = j - 1
        Else If (i .Le. nd-j) Then
           kl = i
           i = i + 1
        Else
           kr = j
           j = j - 1
        End If
     End Do
  End If
  ! adjust kl and kr until the minimal gap ratio is satisfied
  kl0 = kl
  kr0 = kr
  gmin = trl_min_gap_ratio(info, nd, tind, res)
  Do While (kl+minsep.Lt.kr .And. gap_ratio(Max(1,kl),Min(nd,kr)).Lt.gmin)
     If (info%lohi .Lt. 0) Then
        kl = kl + 1
     Else If (info%lohi .Gt. 0) Then
        kr = kr - 1
     Else If (res(kl) .Lt. res(kr)) Then
        kl = kl + 1
     Else
        kr = kr + 1
     End If
  End Do
  ! make sure nearly degenerate Ritz pairs are included
  ! lambda(kl)+r(kl) > lambda(j) and 
  !                lambda(kl)-r(kl) > lambda(j)-r(j)
  If (info%lohi .Gt. 0) Then
     i = kr0-1
     Do While (i.Gt.kl+minsep .And. lambda(kr)-res(kr).Lt.lambda(i) .And.&
          & lambda(kr)+res(kr).Lt.lambda(i)+res(i))
        i = i - 1
     End Do
     kr = Min(kr, i+1)
  Else
     i = kl0 + 1
     Do While (i.Lt.kr-minsep .And. lambda(kl)+res(kl).Gt.lambda(i) .And.&
          & lambda(kl)-res(kl).Gt.lambda(i)-res(i))
        i = i + 1
     End Do
     kl = Max(kl, i-1)
  End If
Contains
  Function gap_ratio(i,j) Result(gamma)
    Integer, Intent(in) :: i, j
    Real(8) :: gamma
    gamma = (lambda(i) - lambda(tind)) / (lambda(j) - lambda(tind))
    Return
  End Function gap_ratio
End Subroutine trl_restart_small_res
!!!
! search throught all pairs of (kl, kr) for the one with the maximum
! gap ratio for the next Ritz pair(target)
!
! This is an optimized version of the original version.  It only search
! through nd number once. (Only single loop!)
!!!
Subroutine trl_restart_max_gap_ratio(nd, tind, kept, lambda, res, info, kl, kr)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T), Intent(in) :: info
  Integer, Intent(in) :: kept, nd, tind
  Integer, Intent(inout) :: kl, kr
  Real(8), Intent(in) :: lambda(nd), res(nd)
  Interface
     Subroutine trl_restart_search_range(nd, lambda, res, info, ncl,&
          & ncr, lohi, tind, klm, krm)
       Use trl_info
       Implicit None
       Type(TRL_INFO_T), Intent(in) :: info
       Integer, Intent(in) :: nd, ncl, ncr, tind
       Integer, Intent(out) :: klm, krm, lohi
       Real(8), Intent(in) :: lambda(nd), res(nd)
     End Subroutine trl_restart_search_range
     Function trl_min_gap_ratio(info, nd, tind, res) Result(gamma)
       Use trl_info
       Implicit None
       Type(TRL_INFO_T), Intent(in) :: info
       Integer, Intent(in) :: nd, tind
       Real(8), Intent(in) :: res(nd)
       Real(8) :: gamma
     End Function trl_min_gap_ratio
  End Interface
  !
  ! local variables
  Integer :: i, j, lohi, klm, krm, igap
  Real(8) :: bnd, tmp, gratio
  ! statement function for computing gap ratio
  gratio(i,j) = (lambda(i)-info%trgt)/(lambda(j)-info%trgt)
  ! determine the search range
  Call trl_restart_search_range(nd, lambda, res, info, kl, kr, lohi,&
       & tind, klm, krm)
  kl = klm
  kr = krm
  igap = Max(Min(nd-info%ned, Nint((krm-klm)*0.4D0)), 2)
  If (igap.Gt.2 .And. info%matvec.Gt.info%maxlan) Then
     If (info%clk_op+info%tick_o.Gt.1.0D1*(info%clk_orth&
          &+info%tick_h+info%clk_res+info%tick_r)) Then
        igap = Max(2, nd-kept-1)
     Else
        bnd = trl_min_gap_ratio(info, nd, tind, res)
        If (info%crat .Lt. bnd) igap = Max(2, nd-kept-1)
     End If
  End If
  If (kl+igap .Gt. kr) Return
  ! two cases depending on lohi
  If (lohi .Gt. 0) Then
     ! target is at the high end of spectrum
     bnd = gratio(kr, kl)
     Do i = klm, krm-igap
        j = i + igap
        tmp = gratio(j,i)
        If (tmp .Gt. bnd) Then
           kl = i
           kr = j
           bnd = tmp
        End If
     End Do
  Else
     bnd = gratio(kl, kr)
     Do i = klm, krm-igap
        j = i + igap
        tmp = gratio(i,j)
        If (tmp .Gt. bnd) Then
           kl = i
           kr = j
           bnd = tmp
        End If
     End Do
  End If
End Subroutine trl_restart_max_gap_ratio
!!!
! search for a pair (kl, kr) such that the reduction in residual norm
! of the target (info%trgt) will be the largest before next restart
! The merit function is [gap-ratio * (m-k)]
!!!
Subroutine trl_restart_max_progress(nd, tind, kept, lambda, res, info, kl, kr)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T), Intent(in) :: info
  Integer, Intent(in) :: kept, nd, tind
  Integer, Intent(inout) :: kl, kr
  Real(8), Intent(in) :: lambda(nd), res(nd)
  Interface
     Subroutine trl_restart_search_range(nd, lambda, res, info, ncl,&
          & ncr, lohi, tind, klm, krm)
       Use trl_info
       Implicit None
       Type(TRL_INFO_T), Intent(in) :: info
       Integer, Intent(in) :: nd, ncl, ncr, tind
       Integer, Intent(out) :: klm, krm, lohi
       Real(8), Intent(in) :: lambda(nd), res(nd)
     End Subroutine trl_restart_search_range
     Function trl_min_gap_ratio(info, nd, tind, res) Result(gamma)
       Use trl_info
       Implicit None
       Type(TRL_INFO_T), Intent(in) :: info
       Integer, Intent(in) :: nd, tind
       Real(8), Intent(in) :: res(nd)
       Real(8) :: gamma
     End Function trl_min_gap_ratio
  End Interface
  !
  ! local variables
  Integer :: i, j, lohi, klm, krm, igap
  Real(8) :: tmp, ss, merit
  ! merit measure the factor of residual norm reduction
  merit(i,j) = (lambda(i)-info%trgt) * Abs(j-i) / (lambda(j)-info%trgt)
  ! determine the search range
  Call trl_restart_search_range(nd, lambda, res, info, kl, kr, lohi,&
       & tind, klm, krm)
  ! perform the brute-force search
  kl = klm
  kr = krm
  igap = Max(Min(nd-info%ned, Nint((krm-klm)*0.4D0)), 2)
  If (igap.Gt.2 .And. igap+kept.Gt.nd .And. info%crat.Gt.0.0D0) Then
     ss = trl_min_gap_ratio(info, nd, tind, res)
     If (ss .Gt. info%crat) igap = Max(2, nd-kept-1)
  End If
  If (lohi .Gt. 0) Then
     ss = merit(kr, kl)
     Do i = klm, krm-igap
        Do j = i+igap, krm
           tmp = merit(j, i)
           If (tmp.Gt.ss) Then
              ss = tmp
              kl = i
              kr = j
           End If
        End Do
     End Do
  Else
     ss = merit(kl, kr)
     Do i = klm, krm-igap
        Do j = i+igap, krm
           tmp = merit(i,j)
           If (tmp.Gt.ss) Then
              ss = tmp
              kl = i
              kr = j
           End If
        End Do
     End Do
  End If
End Subroutine trl_restart_max_progress
!!!
! search for a pair (kl, kr) such that the reduction in residual norm
! of the target (info%trgt) will be the largest before next restart
! the merit function is [sqrt(gap ratio) * (m-k)]
!!!
Subroutine trl_restart_max_reduction(nd, tind, kept, lambda, res, info, kl, kr)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T), Intent(in) :: info
  Integer, Intent(in) :: kept, nd, tind
  Integer, Intent(inout) :: kl, kr
  Real(8), Intent(in) :: lambda(nd), res(nd)
  Interface
     Subroutine trl_restart_search_range(nd, lambda, res, info, ncl,&
          & ncr, lohi, tind, klm, krm)
       Use trl_info
       Implicit None
       Type(TRL_INFO_T), Intent(in) :: info
       Integer, Intent(in) :: nd, ncl, ncr, tind
       Integer, Intent(out) :: klm, krm, lohi
       Real(8), Intent(in) :: lambda(nd), res(nd)
     End Subroutine trl_restart_search_range
     Function trl_min_gap_ratio(info, nd, tind, res) Result(gamma)
       Use trl_info
       Implicit None
       Type(TRL_INFO_T), Intent(in) :: info
       Integer, Intent(in) :: nd, tind
       Real(8), Intent(in) :: res(nd)
       Real(8) :: gamma
     End Function trl_min_gap_ratio
  End Interface
  !
  ! local variables
  Integer :: i, j, lohi, klm, krm, igap
  Real(8) :: tmp, ss, merit
  ! merit measure the factor of residual norm reduction
  merit(i,j) = Sqrt((lambda(i)-info%trgt)/(lambda(j)-info%trgt)) * Abs(j-i)
  ! determine the search range
  Call trl_restart_search_range(nd, lambda, res, info, kl, kr, lohi,&
       & tind, klm, krm)
  ! perform the brute-force search
  kl = klm
  kr = krm
  igap = Max(Min(nd-info%ned, Nint((krm-klm)*0.4D0)), 2)
  If (igap.Gt.2 .And. igap+kept.Gt.nd .And. info%crat.Gt.0.0D0) Then
     ss = trl_min_gap_ratio(info, nd, tind, res)
     If (ss .Gt. info%crat) igap = Max(2, nd-kept-1)
  End If
  If (lohi .Gt. 0) Then
     ss = merit(kr, kl)
     Do i = klm, krm-igap
        Do j = i+igap, krm
           tmp = merit(j, i)
           If (tmp.Gt.ss) Then
              ss = tmp
              kl = i
              kr = j
           End If
        End Do
     End Do
  Else
     ss = merit(kl, kr)
     Do i = klm, krm-igap
        Do j = i+igap, krm
           tmp = merit(i,j)
           If (tmp.Gt.ss) Then
              ss = tmp
              kl = i
              kr = j
           End If
        End Do
     End Do
  End If
End Subroutine trl_restart_max_reduction
!!!
! determine the search range -- used by the schemes that performs brute
! -force search
!!!
Subroutine trl_restart_search_range(nd, lambda, res, info, ncl, ncr,&
     & lohi, tind, klm, krm)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T), Intent(in) :: info
  Integer, Intent(in) :: nd, ncl, ncr, tind
  Integer, Intent(out) :: klm, krm, lohi
  Real(8), Intent(in) :: lambda(nd), res(nd)
  !
  ! local variables
  Integer :: j
  Real(8) :: bnd
  klm = Max(ncl,1)
  krm = Min(ncr,nd)
  bnd = info%tol * info%anrm
  lohi = info%lohi
  ! make sure there is some distance between the boundary and the
  ! target Ritz value
  If (info%lohi .Gt. 0) Then
     ! save high end
     krm = Min(Max(info%maxlan-info%ned, (info%maxlan+info%nec)/2),&
          & krm, tind-1)
     Do While (krm+krm.Ge.ncl+ncr .And. res(krm).Le.bnd)
        krm = krm - 1
     End Do
  Else If (info%lohi .Lt. 0) Then
     ! save lower end
     klm = Max(Min(info%ned, (info%maxlan+info%nec)/2), tind+1, klm)
     Do While (klm+klm.Le.ncl+ncr .And. res(klm).Le.bnd)
        klm = klm + 1
     End Do
  Else
     ! save both ends
     If (tind-klm .Lt. krm-tind) Then
        lohi = -1
        klm = tind + 1
     Else
        lohi = 1
        krm = tind - 1
     End If
     j = info%locked + klm + nd - krm + 1
     If (j .Gt. 0) Then
        j = j / 2
        klm = klm - j
        krm = krm + j
     End If
  End If
End Subroutine trl_restart_search_range
!!!
! try to determine the minimal gap ratio need to compute all wanted
! eigenvalues 
!!!
Function trl_min_gap_ratio(info, nd, tind, res) Result(gamma)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T), Intent(in) :: info
  Integer, Intent(in) :: nd, tind
  Real(8), Intent(in) :: res(nd)
  Real(8) :: gamma
  !
  gamma = info%maxmv * (info%nec + 1.0D0) / info%ned - info%matvec
  If (gamma .Lt. info%klan) Then
     gamma = Max((info%maxmv-info%matvec)/Dble(info%ned-info%nec), 2.0D0)
  End If
  gamma = Min(Log(res(tind)/(info%tol*info%anrm))/gamma, 0.5D0)
  Return
End Function trl_min_gap_ratio