trlcore.f90

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

          title = 'Ritz values saved (may not be ordered)..'
       End If
       Call trl_print_real(info, title, alpha(1:kept+locked))
       title = 'Residual norms of the saved Ritz pairs..'
       betrot(1:kept+locked) = Abs(beta(1:kept+locked))
       Call trl_print_real(info, title, betrot(1:kept+locked))
    End If
    If (info%verbose .Gt. 7) Then
       Do j1 = 1, Min(kept, info%verbose)
          Do j2 = 1, j1
             wrk2(j2) = Dot_product(yy(j2*jml-jml+1:j2*jml),&
                  & yy(j1*jml-jml+1:j1*jml))
          End Do
          wrk2(j1) = wrk2(j1) - 1
          Write(title, *) 'Orthogonality level of y(', j1, ') ..'
          Call trl_print_real(info, title, wrk2(1:j1))
       End Do
    End If
    If (info%verbose .Gt. 10) Then
       Do j1 = 1, Min(kept, info%verbose)
          Write(title,*) 'eigenvector ', j1, ' of Q''AQ..'
          Call trl_print_real(info, title, yy((j1-1)*jml+1:j1*jml))
       End Do
    End If
    If (info%verbose .Gt. 10) Then
       j1n = Min(nrow, info%verbose)
       Do j1 = 1, Min(kept+locked, mev)
          Write(title,*) 'Ritz vector ', j1, ' (1:', j1n, ') ..'
          Call trl_print_real(info, title, evec(1:j1n,j1))
       End Do
    End If
  End Subroutine print_restart_state
  ! print the final state
  Subroutine print_final_state
    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 = 'Final eigenvalues  (in ascending order)..'
    Call trl_print_real(info, title, eval(1:kept))
    If (info%verbose .Gt. 4) Then
       title = 'Final residual norms..'
       Call trl_print_real(info, title, beta(1:kept))
    End If
    If (info%verbose .Gt. 8) Then
       Do j1 = 1, Min(kept, info%verbose)
          Write(title,*) 'Eigenvector ', j1, ' of Q''AQ ..'
          Call trl_print_real(info, title, yy((j1-1)*jml+1:j1*jml))
       End Do
    End If
    If (info%verbose .Gt. 10) Then
       j1n = Min(nrow, info%verbose)
       Do j1 = 1, Min(kept, mev)
          Write(title,*) 'Ritz vector ', j1, ' (1:', j1n, ') ..'
          Call trl_print_real(info, title, evec(1:j1n,j1))
       End Do
    End If
  End Subroutine print_final_state
  ! do check pointing
  Subroutine write_checkpoint
    Interface
       ! write the checkpoint files
       Subroutine trl_write_checkpoint(cp_io, filename, nrow, alpha,&
            & beta, evec, lde, me, base, ldb, nb, ierr)
         Implicit None
         Character(*), Intent(in) :: filename
         Integer, Intent(in) :: cp_io, nrow, ldb, lde, me, nb
         Real(8), Intent(in) :: alpha(me+nb-1), beta(me+nb-1)
         Real(8), Intent(in) :: evec(lde,me), base(ldb,nb)
         Integer, Intent(out) :: ierr
       End Subroutine trl_write_checkpoint
       ! find the minimum of the flag on all PEs
       Function trl_sync_flag(mpicom, inflag) Result(outflag)
         Implicit None
         Integer :: outflag
         Integer, Intent(in) :: mpicom, inflag
       End Function trl_sync_flag
       Function trl_pe_filename(base, my_rank, npe) Result(filename)
         Implicit None
         Integer, Intent(in) :: my_rank, npe
         Character(*), Intent(in) :: base
         Character(132) :: filename
       End Function trl_pe_filename
    End Interface
    Integer :: ii, c1, c2
    title = trl_pe_filename(info%cpfile, info%my_pe, info%npes)
    Call System_clock(c1)
    Call trl_write_checkpoint(info%cpio, title, nrow, alpha,&
         & beta, evec, lde, j1n, base, ldb, j2n, ii)
    Call System_clock(c2)
    If (c2.Gt.c1) Then
       info%clk_out = info%clk_out + (c2 - c1)
    Else
       info%clk_out = info%clk_out + ((info%clk_max-c1)+c2)
    End If
    info%wrds_out = info%wrds_out + jnd*(nrow+nrow+2) + nrow + 2
    info%stat = trl_sync_flag(info%mpicom, ii)
  End Subroutine write_checkpoint
  ! dump important variables when return due to error
  Subroutine log_error_state
    Interface
       ! print real array for debugging
       Subroutine trl_print_real(info, title, array)
         Use trl_info
         Type(TRL_INFO_T), Intent(in) :: info
         Character(*), Intent(in) :: title
         Real(8), Dimension(:), Intent(in) :: array
       End Subroutine trl_print_real
       ! print integer array for debugging
       Subroutine trl_print_int(info, title, array)
         Use trl_info
         Type(TRL_INFO_T), Intent(in) :: info
         Character(*), Intent(in) :: title
         Integer, Dimension(:), Intent(in) :: array
       End Subroutine trl_print_int
       ! print a time stamp
       Subroutine trl_time_stamp(iou)
         Implicit None
         Integer, Intent(in) :: iou
       End Subroutine trl_time_stamp
       ! a short of the print info function that can be called by
       ! indivadual PE
       Subroutine trl_terse_info(info, iou)
         Use trl_info
         Implicit None
         Type(TRL_INFO_T), Intent(in) :: info
         Integer, Intent(in) :: iou
       End Subroutine trl_terse_info
    End Interface
    Call trl_time_stamp(info%log_io)
    title = 'Dumping the content of the variables on error..'
    iwrk(1) = info%stat
    Call trl_print_int(info, title, iwrk(1:1))
    i2 = info%log_io
    Call trl_terse_info(info, i2)
    Write (i2,*) 'This Lanczos iteration started With ', kept,&
         & ' vectors.'
    Write (i2,*) 'There are ', jnd, '(', j1, ', ', j2,&
         & ') Lanczos vectors currently.'
    If (jnd .Ne. j1+j2) jnd = j1 + j2
    If (jnd.Lt.0 .Or. jnd.Gt.info%klan) jnd = 0
    title = 'Content of eval ..'
    Call trl_print_real(info, title, eval)
    If (jnd .Gt. 0) Then
       Write (title, *) 'Alpha(1:', jnd, ') .. '
       Call trl_print_real(info, title, alpha(1:jnd))
       title(1:5) = ' Beta'
       Call trl_print_real(info, title, beta(1:jnd))
       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 If
    If (j1 .Gt. 0) Then
       title = 'the First row of evec ..'
       Call trl_print_real(info, title, evec(1, 1:j1))
       Write (title, *) 'row ', nrow, ' of evec ..'
       Call trl_print_real(info, title, evec(nrow, 1:j1))
    End If
    If (j2 .Gt. 0) Then
       title = 'the First row of base ..'
       Call trl_print_real(info, title, base(1,1:j2))
       Write (title, *) 'row ', nrow, ' of base ..'
       Call trl_print_real(info, title, base(nrow,1:j2))
    End If
    If (Associated(qb)) Then
       Write (title, *) 'Content of qb (q_', jnd-1, ') ..'
       Call trl_print_real(info, title, qb)
    End If
    If (Associated(qa)) Then
       Write (title, *) 'Content of qa (q_', jnd, ') ..'
       Call trl_print_real(info, title, qa)
    End If
    If (Associated(rr)) Then
       Write (title, *) 'Content of rr (residual == q_', jnd+1, ') ..'
       Call trl_print_real(info, title, rr)
    End If
    If (info%my_pe.Eq.0 .And. info%log_io.Ne.6) Then
       Write (*, *) 'TRLanczos returned with error ', info%stat
       If (info%verbose .Gt. 0) Then
          Write (*, *) 'Contents of most variables are dumped to ',&
               & 'log file ', Trim(info%log_file)
       Else
          Write (*, *) 'Contents of most variables are dumped to ',&
               &'IO unit # ', info%log_io
       End If
    End If
  End Subroutine log_error_state
End Subroutine trlanczos
!!!
! check to make sure the initial guess vector contains valid nonzero
! numbers if not fill with random numbers
! this routine will also read the checkpoint files to restore the
! previous state of the Lancozs iterations
!!!
Subroutine trl_initial_guess(nrow, evec, lde, mev, base, ldb, nbas,&
     & alpha, beta, j1, j2, info, wrk, lwrk)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T) :: info
  Integer, Intent(in) :: nrow, lde, ldb, mev, nbas, lwrk
  Integer, Intent(out) :: j1, j2
  Real(8) :: wrk(lwrk), base(ldb,nbas), evec(lde,mev)
  Real(8) :: alpha(mev+nbas-1), beta(mev+nbas-1)
  Interface
     ! read the checkpointfile
     Subroutine trl_read_checkpoint(cp_io, filename, nrow, evec, lde,&
          & mev, j1, base, ldb, nbas, j2, alpha, beta, ierr)
       Implicit None
       Character(*), Intent(in) :: filename
       Integer, Intent(in) :: cp_io, nrow, lde, mev, ldb, nbas
       Integer, Intent(out) :: j1, j2, ierr
       Real(8), Intent(out) :: alpha(mev+nbas-1), beta(mev+nbas-1)
       Real(8), Intent(out) :: base(ldb,nbas), evec(lde,mev)
     End Subroutine trl_read_checkpoint
     ! the orthogonalization routine -- Classical Gram-Schmidt
     Subroutine trl_cgs(mpicom, myid, nrow, v1, ld1, m1, v2, ld2,&
          & m2, rr, rnrm, alpha, north, nrand, wrk, ierr)
       Integer, Intent(inout) :: north, nrand, ierr
       Integer, Intent(in) :: mpicom, myid, nrow, ld1, ld2, m1, m2
       Real(8), Intent(inout) :: rnrm, alpha
       Real(8), Intent(in) :: v1(ld1,m1), v2(ld2*m2)
       Real(8), Intent(inout) :: rr(nrow), wrk(m1+m2+m1+m2)
     End Subroutine trl_cgs
     ! generate filenames for each PE
     Function trl_pe_filename(base, my_rank, npe) Result(filename)
       Implicit None
       Integer, Intent(in) :: my_rank, npe
       Character(*), Intent(in) :: base
       Character(132) :: filename
     End Function trl_pe_filename
     ! find the minimum of the flag on all PEs
     Function trl_sync_flag(mpicom, inflag) Result(outflag)
       Implicit None
       Integer :: outflag
       Integer, Intent(in) :: mpicom, inflag
     End Function trl_sync_flag
     Subroutine trl_g_sum(mpicom, n, xx, wrk)
       Implicit None
       Integer, Intent(in) :: mpicom, n
       Real(8), Intent(out) :: wrk(n)
       Real(8), Intent(inout) :: xx(n)
     End Subroutine trl_g_sum
     Subroutine trl_check_orth(info, nrow, v1, ld1, j1, v2, ld2, j2, wrk, lwrk)
       Use trl_info
       Implicit None
       Integer, Intent(in) :: nrow, ld1, ld2, j1, j2, lwrk
       Type(TRL_INFO_T), Intent(in) :: info
       Real(8), Intent(in) :: v1(ld1,j1), v2(ld2,j2)
       Real(8) :: wrk(lwrk)
     End Subroutine trl_check_orth
  End Interface
  !
  ! local variable
  !
  Integer :: i, ii, j, nran, north
  Real(8) :: tmp, rnrm
  Character(132) :: file
  Integer, Dimension(:), Allocatable :: rseed
  !
  ! generate random seeds based on current clock ticks
  Call Random_seed(SIZE=ii)    ! find out the size of seed
  Allocate(rseed(ii))
  Call Random_seed(GET=rseed)  ! get the current seeds
  i = Mod(info%my_pe, ii) + 1
  j = Max(Maxval(rseed), info%nloc)
  Call System_clock(rseed(i))  ! get the current clock ticks
  ! just in case the clocks are absolutely synchronized
  If (info%my_pe .Gt. 0) Then
     rseed(i) = rseed(i) - Sign(Int(info%my_pe*Sqrt(Dble(j))), rseed(i))
  End If
  Call Random_seed(PUT=rseed)  ! initialize the random number generator
  Deallocate(rseed)
  !
  j = info%nec+1
  If (info%guess .Gt. 1) Then
     ! retrieve a check-point file
     i = info%cpio
     If (info%oldcpf .Ne. '') Then
        file = trl_pe_filename(info%oldcpf, info%my_pe, info%npes)
     Else
        file = trl_pe_filename(info%cpfile, info%my_pe, info%npes)
     End If
     Call System_clock(ii)
     Call trl_read_checkpoint(info%cpio, file, nrow, evec(1,j), lde,&
          & mev-info%nec, j1, base, ldb, nbas, j2, alpha(j), beta(j),&
          & i)
     info%stat = trl_sync_flag(info%mpicom, i)
     Call System_clock(i)
     If (i.Gt.ii) Then
        info%clk_in = i - ii
     Else
        info%clk_in = (info%clk_max - ii) + i
     End If
     info%wrds_in = (j1+j2)*(nrow+nrow+2) + nrow + 2
     j1 = j1 + info%nec
     If (info%stat .Ne. 0) Return
  Else
     If (info%guess .Le. 0) Then
        ! generate an arbitrary initial starting vector
        ! if (info%guess == 0), use the vector [1, 1, ...]^T
        ! else perturb some random locations of the above vector
        evec(1:nrow,j) = 1.0D0
        nran = Min(1-info%guess, lwrk)
        nran = 2*(nran/2)
        If (nran.Gt.0 .And. nran.Lt.nrow) Then
           Call Random_number(wrk(1:nran))
           Do i = 1, nran-1, 2
              ii = Int(nrow*wrk(i))+1
              evec(ii,j) = evec(ii,j) + wrk(i+1) - 0.5D0
           End Do
           info%flop = info%flop + nran + nran
        Else If (nran .Ge. nrow) Then
           Call random_vector
        End If
     End If
     j1 = info%nec
     j2 = 0
  End If
  tmp = 0.0
  ! make sure the norm of the next vector can be computed
  wrk(1) = Dot_product(evec(1:nrow,j), evec(1:nrow,j))
  Call trl_g_sum(info%mpicom, 1, wrk(1), wrk(2))
  info%flop = info%flop + nrow+nrow
  If (wrk(1).Ge.Tiny(tmp) .And. wrk(1).Le.Huge(tmp)) Then
     ! set rnrm to let trl_CGS normalize evec(1:nrow, j)
     rnrm = Sqrt(wrk(1))
  Else
     Call random_vector
  End If
  !
  ! orthogonalize initial guess against all existing vectors
  !
  i = 0
  tmp = 1.0D0
  nran = info%nrand
  north = info%north
  If (j1 .Lt. mev) Then
     Call trl_cgs(info%mpicom, info%my_pe, nrow, evec, lde, j1,&
          & base, ldb, 0, evec(1,j1+1), rnrm, tmp, i,&
          & info%nrand, wrk, info%stat)
  Else If (j2 .Le. 0) Then
     Call trl_cgs(info%mpicom, info%my_pe, nrow, evec, lde, j1,&
          & evec, lde, 0, base, rnrm, tmp, i, info%nrand, wrk,&
          & info%stat)
  Else
     Call trl_cgs(info%mpicom, info%my_pe, nrow, evec, lde, j1,&
          & base, ldb, j2, base(1,j2+1), rnrm, tmp, i, info%nrand,&
          & wrk, info%stat)
  End If
  info%flop = info%flop + 4*nrow*((info%north-north)*(j1+j2) + &
       & info%nrand-nran) + nrow
  If (info%verbose.Gt.6) Then
     If (j1 .Lt. mev) Then
        i = j1 + 1
        ii = j2
     Else
        i = j1
        ii = j2 + 1
     End If
     Call trl_check_orth(info, nrow, evec, lde, j1, base, ldb, ii, wrk,&
          & lwrk)
  End If
  Return
Contains
  ! an internal subroutine to generate random vector
  Subroutine random_vector
    Call Random_number(evec(1:nrow, j))
    evec(1:nrow, info%nec+1) = evec(1:nrow, info%nec+1) +&
         & evec(1:nrow, info%nec+1) +&
         & Cshift(evec(1:nrow, info%nec+1), 1) +&
         & Cshift(evec(1:nrow, info%nec+1), -1)
    info%nrand = info%nrand + 1
    info%flop = info%flop + 4*nrow
  End Subroutine random_vector
End Subroutine trl_initial_guess
!!!
! full re-orthogonalization for TRLANCZOS
!
!\Algorithm
!