trlcore.f90

用于计算矩阵的特征值,及矩阵的其他运算.可以用与稀疏矩阵

        info%matvec = info%matvec+1
        ! compute alpha(jnd)
        alpha(jnd) = Dot_product(qa, rr)
        ! the Lanczos orthogonalization
        Call trl_g_sum(info%mpicom, 1, alpha(jnd:jnd), wrk2)
        rr = rr - alpha(jnd)*qa - beta(jnd-1)*qb
        info%flop = info%flop + 6*nrow
        ! re-orthogonalization, compute beta(jnd)
        info%flop_h = info%flop_h - info%flop
        Call System_clock(clk1)
        Call trl_orth(nrow, evec, lde, j1, base, ldb, j2, rr, kept,&
             & alpha, beta, wrk2, lwrk2, info)
        If (info%verbose.Gt.8 .Or. info%stat.Ne.0) Then
           ! check orthogonality after the initilization step
           Call trl_check_orth(info, nrow, evec, lde, j1n, base, ldb, j2n,&
                & wrk2, lwrk2)
        End If
        Call add_clock_ticks(info%clk_orth, info%tick_h)
        info%flop_h = info%flop_h + info%flop
        If (info%stat.Ne.0) Goto 888
        alfrot(jnd) = alpha(jnd)
        betrot(jnd) = beta(jnd)
        If (info%verbose .Gt. 4) Call print_alpha_beta
        !
        ! perform convergence test once in a while
        If (info%matvec.Ge.next_test .And. (beta(jnd).Ne.zero .Or. &
             & info%matvec.Gt.info%maxlan)) Then
           If (info%verbose .Gt. 5) Call print_all_alpha_beta
           lambda => wrk2(1:jnd)
           res => wrk2(jnd+1:jnd+jnd)
           Call trl_get_eval(jnd, locked, alfrot, betrot, lambda, res,&
                & wrk2(jnd+jnd+1:lwrk2), lwrk2-jnd-jnd, info%stat)
           If (info%stat .Ne. 0) Goto 888
           If (info%verbose .Gt. 2) Call print_lambda_res
           i1 = Min(mev, jnd)
           eval(1:i1) = wrk2(1:i1)
           Call trl_convergence_test(jnd, lambda, res, info,&
                & wrk2(jnd+jnd+1:4*jnd))
           ! decide when to perform the next test
           If (info%nec .Lt. info%ned .And. info%nec.Gt.0) Then
              next_test = Dble(info%ned * info%matvec) / Dble(info%nec)
           Else If (info%nec.Eq.0) Then
              next_test = next_test + next_test
              If (info%maxlan .Eq. info%ntot) Then
                 next_test = Ceiling(0.5*(info%maxlan + info%matvec))
              End If
           End If
           If (info%verbose .Gt. 0) Call trl_print_progress(info)
        End If
     !!**********************************************************!!
     End Do ! inner (regular Lanczos three-term recurrence) loop !!
     !!**********************************************************!!
     !
     ! error checking for debugging use
     !
     lambda => wrk2(1:jnd)
     res => wrk2(jnd+1:jnd+jnd)
     If (info%verbose.Gt.6) Then
        wrk2 => wrk2(jnd+jnd+1:lwrk2)
        i2 = lwrk2 - jnd - jnd
        Call trl_check_orth(info, nrow, evec, lde, j1n, base, ldb, j2n,&
             & wrk2, i2)
        If (info%verbose.Gt.7) Then
           Call trl_check_recurrence(op, info, nrow, evec, lde, j1n,&
                & base, ldb, j2n, kept, alpha, beta, wrk2, i2)
        End If
     End If
     !
     ! convert the integer counters to floating-point counters
     !
     i2 = info%clk_max / 4
     If (info%flop .Gt. i2) Then
        info%rflp = info%rflp + info%flop
        info%flop = 0
     End If
     If (info%flop_h .Gt. i2) Then
        info%rflp_h = info%rflp_h + info%flop_h
        info%flop_h = 0
     End If
     If (info%flop_r .Gt. i2) Then
        info%rflp_r = info%rflp_r + info%flop_r
        info%flop_r = 0
     End If
     If (info%clk_op .Gt. i2) Then
        info%tick_o = info%tick_o + info%clk_op
        info%clk_op = 0
     End If
     If (info%clk_orth .Gt. i2) Then
        info%tick_h = info%tick_h + info%clk_orth
        info%clk_orth = 0
     End If
     If (info%clk_res .Gt. i2) Then
        info%tick_r = info%tick_r + info%clk_res
        info%clk_res = 0
     End If
     info%flop_r = info%flop_r - info%flop
     !
     ! *** Determine whether to restart ***
     ! compute the Ritz values and Ritz vectors if they are not up to
     ! date
     !
     Call System_clock(clk1)
     prek = kept
     jml = jnd - locked
     i2 = kept - locked + 1
     If (info%nec .Lt. info%ned) Then
        ! need to compute the updated Ritz values and residual norms
        wrk2 => wrk(4*info%maxlan+2*jnd+1:lwrk-i2*i2)
        lwrk2 = lwrk-i2*i2-4*info%maxlan-2*jnd
        If (lwrk2 .Lt. 3*jnd) Then
           info%stat = -12
           Goto 888
        End If
        If (info%verbose .Gt. 5) Call print_all_alpha_beta
        Call trl_get_eval(jnd, locked, alfrot, betrot,&
             & lambda, res, wrk2, lwrk2, info%stat)
        If (info%stat .Ne. 0) Goto 888
        If (info%verbose .Gt. 2) Call print_lambda_res
        Call trl_convergence_test(jnd, lambda, res, info, wrk2)
        If (info%verbose .Gt. 0) Call trl_print_progress(info)
     End If
     !
     ! compute the Ritz vectors
     ! decide how many vectors to save if restart
     If (info%nec.Lt.info%ned .And. info%matvec.Lt.info%maxmv) Then
        ! prepare to restart
        Call trl_shuffle_eig(jml, info%klan-locked, lambda(locked&
             &+1:locked+jml), res(locked+1:locked+jml), info, kept)
        ! compute eigenvectors using dstein (inverse interations)
        If (kept*3 .Lt. jml) Then
           i1 = 4*info%maxlan+kept*jml+jnd
           yy => wrk(4*info%maxlan+jnd+1:i1)
           wrk2 => wrk(i1+1:lwrk-i2*i2)
           lwrk2 = lwrk - i1 - i2*i2
           Call trl_get_tvec(jml, alfrot(locked+1:locked+jml),&
                & betrot(locked+1:locked+jml), 0, i2, rot, kept,&
                & lambda(locked+1:locked+jml), yy, iwrk, wrk2, lwrk2,&
                & info%stat)
           If (info%stat.Eq.0)&
                & alpha(1:locked+kept) = lambda(1:locked+kept)
        End If
        ! compute eigenvectors using dsyev (QR)
        If (kept*3.Ge.jml .Or. info%stat.Ne.0) Then
!0918        If (info%stat.Ne.0) Then
           alfrot(1:locked+kept) = lambda(1:locked+kept)
           i1 = 4*info%maxlan+jml*jml
           yy => wrk(4*info%maxlan+1:i1)
           wrk2 => wrk(i1+1:lwrk)
           lwrk2 = lwrk - i1
           Call trl_get_tvec_a(jml, prek-locked, alpha(locked+1:locked&
                &+jml), beta(locked+1:locked+jml), kept, alfrot(locked&
                &+1:locked+jml), yy, wrk2, lwrk2, iwrk, info%stat)
        End If
        If (info%stat .Ne. 0) Goto 888
        beta(locked+1:locked+kept) = yy(jml:jml*kept:jml)*betrot(jnd)
        If (jml.Gt.info%ned+(info%ned/5+6)) Then
           Call trl_set_locking(jml, kept, alpha(locked+1:locked+kept),&
                & beta(locked+1:locked+kept), yy, info%anrm, i2)
        Else
           i2 = 0
        End If
        !
        ! generate Ritz vectos, reclaim the space pointed by ROT
        i1 = 4*info%maxlan + kept*jml + jnd
        wrk2 => wrk(i1+1:lwrk)
        lwrk2 = lwrk - i1
        Call trl_ritz_vectors(nrow, locked, kept, yy, jml,&
             & evec, lde, j1, base, ldb, j2, wrk2, lwrk2)
        info%flop = info%flop + 2*nrow*jml*kept
        If (info%verbose .Gt. 0) Call print_restart_state
        ! reset the counters and indices to the correct values for
        ! restarting
        kept = kept + locked
        locked = locked + i2
        info%locked = locked
        jnd = kept
        If (jnd .Le. mev) Then
           j1 = jnd
           j2 = 0
        Else
           j1 = mev
           j2 = jnd - mev
        End If
        If (info%nec .Gt. 0) Then
           next_test = Int(Dble(info%matvec*info%ned)/Dble(info%nec))
        Else
           next_test = next_test + info%maxlan
        End If
        i1=Min(mev, jnd)
        eval(1:i1) = lambda(1:i1)
        If (jnd .Lt. mev) Then
           j1n = j1 + 1
           j2n = 0
           evec(1:nrow, j1n) = rr
        Else
           j1n = mev
           j2n = j2 + 1
           base(1:nrow, j2n) = rr
        End If
        ! write checkpoint files
        If (info%matvec .Ge. chkpnt) Then
           Call write_checkpoint
           chkpnt = chkpnt + info%maxmv/info%cpflag
        End If
     Else
        !
        ! all wanted eigenpairs converged or maximum MATVEC used
        ! sort the eigenvalues in final output order
        kept = Min(info%nec, Max(info%ned,mev-1))
        info%nec = kept
        If (kept.Eq.0) kept = Min(mev-1, info%ned)
        Call trl_sort_eig(jnd, info%lohi, kept, lambda, res)
        eval(1:kept) = lambda(1:kept)
        If (kept*3 .Lt. jnd) Then
           ! eigenvectors of the projection matrix (try inverse
           ! interations)
           i1 = kept*jnd + 4*info%maxlan
           yy => wrk(4*info%maxlan+1:i1)
           wrk2 => wrk(i1+1:lwrk-i2*i2)
           lwrk2 = lwrk - i1 - i2*i2
           Call trl_get_tvec(jnd, alfrot, betrot, locked, i2, rot,&
                & kept, eval, yy, iwrk, wrk2, lwrk2, info%stat)
        End If
        If (kept*3.Ge.jnd .Or. info%stat .Ne. 0) Then
           ! too many eigenvectors or inverse iterations have failed,
           ! try QR
           i1 = 4*info%maxlan+jnd*jnd
           yy => wrk(4*info%maxlan+1:i1)
           wrk2 => wrk(i1+1:lwrk)
           lwrk2 = lwrk - i1
           Call trl_get_tvec_a(jnd, prek, alpha, beta, kept, eval, yy,&
                & wrk2, lwrk2, iwrk, info%stat)
           If (info%stat .Ne. 0) Goto 888
        End If
        alpha(1:kept) = eval(1:kept)
        beta(1:kept) = betrot(jnd)*yy(jnd:kept*jnd:jnd)
        !
        ! generate eigenvectos, reclaim the space pointed by ROT
        i1 = kept*jnd + 4*info%maxlan
        wrk2 => wrk(i1+1:lwrk)
        lwrk2 = lwrk - i1
        Call trl_ritz_vectors(nrow, 0, kept, yy, jnd,&
             & evec, lde, j1, base, ldb, j2, wrk2, lwrk2)
        info%flop = info%flop + 2*nrow*jml*kept
        If (info%verbose .Gt. 1) Call print_final_state
        ! reset the counters and indices to be used by check_orth and
        ! check_recurrence
        jnd = kept
        j1 = kept
        j2 = 0
        If (kept .Lt. mev) Then
        End If
        If (j1 .Lt. mev) Then
           j1n = j1 + 1
           j2n = 0
           evec(1:nrow, j1n) = rr
        Else
           j1n = mev
           j2n = 1
           base(1:nrow, 1) = rr
        End If
        ! write checkpoint files
        If (info%cpflag .Gt. 0) Call write_checkpoint
     End If
     !
     ! check the orthogonality of the basis vectors before restarting
     !
     If (info%verbose.Gt.6) Then
        Call trl_check_orth(info, nrow, evec, lde, j1n, base, ldb, j2n,&
             & wrk2, lwrk2)
        If (info%verbose.Gt.7) Then
           Call trl_check_recurrence(op, info, nrow, evec, lde, j1n,&
                & base, ldb, j2n, kept, alpha, beta, wrk2, lwrk2)
        End If
     End If
     Call add_clock_ticks(info%clk_res, info%tick_r)
     info%flop_r = info%flop_r + info%flop
     info%nloop = info%nloop + 1
  End Do
  !!*********************!!
  !! end of restart_loop !!
  !!*********************!!
  ! write the estimated residual norms to the beginning of WRK
  wrk(1:j1) = Abs(beta(1:j1))
888 If (info%stat.Lt.0 .And. (info%verbose.Gt.0.Or.info%my_pe.Eq.0))&
       & Call log_error_state
  Deallocate(iwrk)
  Return
  ! internal subroutines for printing etc
Contains
  ! add clock ticks to a integer variable, if there is potential for
  ! overflow convert it to floating-point numbers
  Subroutine add_clock_ticks(time, rtime)
    Integer :: time, clk2
    Real(8) :: rtime
    Call System_clock(clk2)
    If (clk2 .Ge. clk1) Then
       clk2 = clk2 - clk1
    Else
       clk2 = clk2 + (info%clk_max - clk1)
    End If
    If (clk2+time .Ge. time) Then
       time = clk2 + time
    Else
       rtime = (rtime + time) + clk2
       time = 0
    End If
  End Subroutine add_clock_ticks
  ! print the current alpha(i) and beta(i)
  Subroutine print_alpha_beta
    Interface
       ! print real array for debugging
       Subroutine trl_print_real(info, title, array)
         Use trl_info
         Implicit None
         Type(TRL_INFO_T), Intent(in) :: info
         Character(*), Intent(in) :: title
         Real(8), Dimension(:), Intent(in) :: array
       End Subroutine trl_print_real
    End Interface
    Write(title, *) 'alpha(',jnd,') ='
    Call trl_print_real(info, title, alpha(jnd:jnd))
    Write(title, *) ' beta(',jnd,') ='
    Call trl_print_real(info, title, beta(jnd:jnd))
  End Subroutine print_alpha_beta
  ! print all computed alpha and beta
  Subroutine print_all_alpha_beta
    Interface
       ! print real array for debugging
       Subroutine trl_print_real(info, title, array)
         Use trl_info
         Implicit None
         Type(TRL_INFO_T), Intent(in) :: info
         Character(*), Intent(in) :: title
         Real(8), Dimension(:), Intent(in) :: array
       End Subroutine trl_print_real
    End Interface
    Write(title, *) 'alfrot(1:', jnd, ')..'
    Call trl_print_real(info, title, alfrot(1:jnd))
    title(1:3) = 'bet'
    Call trl_print_real(info, title, betrot(1:jnd))
  End Subroutine print_all_alpha_beta
  ! print lambda and residual norm
  Subroutine print_lambda_res
    Interface
       ! print real array for debugging
       Subroutine trl_print_real(info, title, array)
         Use trl_info
         Implicit None
         Type(TRL_INFO_T), Intent(in) :: info
         Character(*), Intent(in) :: title
         Real(8), Dimension(:), Intent(in) :: array
       End Subroutine trl_print_real
    End Interface
    title = 'Current eigenvalues..'
    Call trl_print_real(info, title, lambda)
    title = 'Current residual norms..'
    Call trl_print_real(info, title, res)
  End Subroutine print_lambda_res
  ! print current solution status
  Subroutine print_restart_state
    Interface
       ! print integer array for debugging
       Subroutine trl_print_int(info, title, array)
         Use trl_info
         Implicit None
         Type(TRL_INFO_T), Intent(in) :: info
         Character(*), Intent(in) :: title
         Integer, Dimension(:), Intent(in) :: array
       End Subroutine trl_print_int
       ! print real array for debugging
       Subroutine trl_print_real(info, title, array)
         Use trl_info
         Implicit None
         Type(TRL_INFO_T), Intent(in) :: info
         Character(*), Intent(in) :: title
         Real(8), Dimension(:), Intent(in) :: array
       End Subroutine trl_print_real
    End Interface
    iwrk(1) = kept+locked
    iwrk(2) = locked+i2
    title = 'Number of saved and locked Ritz pairs..'
    Call trl_print_int(info, title, iwrk(1:2))
    If (info%verbose.Gt.2) Then
       If (iwrk(2) .Eq. 0) Then
          title = 'Ritz values saved (ascending ordered)..'
       Else