trlan.f90

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

!           Each processor will generate its own debug file with the
!           file name formed by appending the processor number to
!           the string FILE
! The actual name is composed by the routine TRL_PE_FILENAME which is
! in trlaux.f90
!!!
Subroutine trl_set_debug(info, msglvl, iou, file)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T) :: info
  Integer, Intent(in) :: msglvl, iou
  Character(*), Optional :: file
  info%verbose = msglvl
  info%log_io = iou
  If (Present(file)) Then
     info%log_file = file
  End If
End Subroutine trl_set_debug
!!!
! set up the information related to check-pointing
!!!
Subroutine trl_set_checkpoint(info, cpflag, cpio, file)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T) :: info
  Integer, Intent(in) :: cpflag, cpio
  Character(*), Optional :: file
  !
  ! copy the flags
  info%cpflag = cpflag
  info%cpio = cpio
  If (Present(file)) Then
     info%cpfile = file
  End If
End Subroutine trl_set_checkpoint
!!!
! set up parameters related to initial guesses of the Lanczos iterations
!
! It sets the number of eigenvector already converged (initially
! assumed to be zero) and whether the user has provided initial guess
! vector (initially assumed no).  It can also tell TRLan to read
! checkpoint file that is different from the default name.
!
! It uses info%cpio to open the checkpoint files.  If the default value
! of info%cpio is in use, make sure that the routine is called after
! trl_set_checkpoint function is called.
!
! If oldcpf is not specified, TRLAN will use the string info%cpfile as
! the name of the checkpoint file to use.  (INFO%cpfile is the name of
! the output checkpoint files.)
!!!
Subroutine trl_set_iguess(info, nec, iguess, oldcpf)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T) :: info
  Integer, Intent(in) :: iguess, nec
  Character(*), Intent(in), Optional :: oldcpf
  Interface
     ! function to synchronize the flags at output
     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
  Character(132) :: cpf
  ! assign nec and iguess flags to info
  info%nec = nec
  info%guess = iguess
  If (Present(oldcpf)) Then
     info%oldcpf = oldcpf
  Else
     info%oldcpf = ''
  End If
  If (info%oldcpf .Ne. '' .And. info%guess.Gt.1) Then
     ! check to make sure the files exist
     cpf = trl_pe_filename(info%oldcpf, info%my_pe, info%npes)
     Open(info%cpio, file=cpf, status='OLD', form='UNFORMATTED',&
          & iostat=info%stat)
     If (info%stat .Eq. 0) Then
        Close(info%cpio, iostat=info%stat)
        If (info%stat.Ne.0) info%stat = -9
     Else
        info%stat = -8
     End If
     info%stat = trl_sync_flag(info%mpicom, info%stat)
  Else
     info%stat = 0
  End If
End Subroutine trl_set_iguess
!!!
! next function provides an uniform way of printing information
! stored in TRL_INFO_T.  It needs to be called by all PEs.
!
! Arguments
! info  -- the TRL_INFO_T variable to be printed
! mvflop-- the number of floating-operations per MATVEC per PE.  This
!          information has to be supplied by user, otherwise related
!          entries are left blank in the print out.
!
! In parallel environment, when writing to standard outputd device, only
! PE0 will actually write its local summary information.
!
! *** NOTE *** MUST be called on all PEs at the same time ***
!!!
Subroutine trl_print_info(info, mvflop)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T), Intent(in) :: info
  Integer, Intent(in), Optional :: mvflop
  Interface
     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
     Subroutine trl_time_stamp(iou)
       Implicit None
       Integer, Intent(in) :: iou
     End Subroutine trl_time_stamp
     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
  End Interface
  Real(8), Parameter :: zero=0.0D0, one=1.0D0
  Integer :: rate, iot, ierr
  Real(8) :: t_tot, t_op, t_orth, t_res, t_in, t_out, rinv
  Real(8) :: r_tot, r_op, r_orth, r_res, r_in, r_out
  Real(8) :: tmp1(12), tmp2(12)
  Character(132) :: filename
  If (info%clk_rate.Gt.0) Then
     rinv = one / Dble(info%clk_rate)
  Else
     ! get clock rate
     Call System_clock(count_rate=rate)
     rinv = one / Dble(rate)
  Endif
  iot = info%log_io
  If (iot.Le.0 .Or. info%verbose.Le.0) iot = 6
  t_op  = (info%tick_o+info%clk_op) * rinv
  t_tot = (info%tick_t+info%clk_tot) * rinv
  t_res  = (info%tick_r+info%clk_res) * rinv
  t_orth = (info%tick_h+info%clk_orth) * rinv
  t_in = info%clk_in * rinv
  t_out = info%clk_out * rinv
  If (t_op.Ne.zero .And. Present(mvflop)) Then
     If (mvflop .Gt. 0) Then
        r_op = mvflop * info%matvec
     Else
        r_op = zero
     End If
  Else
     r_op = zero
  End If
  If (t_orth .Ne. zero) Then
     r_orth = info%rflp_h + info%flop_h
  Else
     r_orth = zero
  End If
  If (t_res .Ne. zero) Then
     r_res = info%rflp_r+info%flop_r
  Else
     r_res = zero
  End If
  If (r_op .Gt. zero) Then
     r_tot = Dble(mvflop)*Dble(info%matvec)+info%rflp+&
          &Dble(info%flop)
  Else
     r_tot = zero
  End If
  If (info%clk_in .Gt. 0) Then
     r_in = 8D0 * info%wrds_in
  Else
     r_in = zero
  End If
  If (info%clk_out .Gt. 0) Then
     r_out = 8D0 * info%wrds_out
  Else
     r_out = zero
  End If
  tmp1(1) = t_tot
  tmp1(2) = t_op
  tmp1(3) = t_orth
  tmp1(4) = t_res
  tmp1(5) = t_in
  tmp1(6) = t_out
  tmp1(7) = r_tot
  tmp1(8) = r_op
  tmp1(9) = r_orth
  tmp1(10) = r_res
  tmp1(11) = r_in
  tmp1(12) = r_out
  Call trl_g_sum(info%mpicom, 12, tmp1, tmp2)
  If (iot.Eq.info%log_io .And. iot.Ne.6) Then
     filename = trl_pe_filename(info%log_file, info%my_pe, info%npes)
     Open(iot, file=filename, status='UNKNOWN', position='APPEND',&
          & action='WRITE', iostat=ierr)
     If (ierr .Ne. 0) iot = 6
  End If
  If (iot.Eq.6 .And. info%my_pe.Gt.0) Return
  Call trl_time_stamp(iot)
  If (info%npes .Gt. 1) Then
     Write (iot, *) 'TRLAN execution summary (exit status =', info%stat,&
          & ') on PE ', info%my_pe
  Else
     Write (iot, *) 'TRLAN execution summary (exit status =', info%stat,&
          & ')'
  End If
  If (info%lohi.Gt.0) Then
     Write (iot, FMT=100) 'LARGEST', info%nec, info%ned
  Else If (info%lohi.Lt.0) Then
     Write (iot, FMT=100) 'SMALLEST', info%nec, info%ned
  Else
     Write (iot, FMT=100) 'EXTREME', info%nec, info%ned
  End If
  Write (iot, FMT=200) 'Times the operator is applied:', info%matvec,&
       & info%maxmv
  Write (iot, FMT=300) info%nloc, info%my_pe, info%ntot
  Write (iot, FMT=400) 'Convergence tolerance:', info%tol,&
       & info%tol*info%anrm
  Write (iot, FMT=500) 'Maximum basis size:', info%maxlan
  Write (iot, FMT=500) 'Restarting scheme:', info%restart
  Write (iot, FMT=500) 'Number of re-orthogonalizations:', info%north
  Write (iot, FMT=500) 'Number of (re)start loops:', info%nloop
  If (info%nrand .Gt. 0) Then
     Write (iot, FMT=500) 'Number of random vectors used:', info%nrand
  Endif
  If (info%npes .Gt. 1) Then
     Write (iot, FMT=500) 'Number of MPI processes:', info%npes
  Endif
  Write (iot, FMT=500) 'Number of eigenpairs locked:', info%locked
  If (t_op .Gt. zero) Then
     Write (iot, FMT=610) 'OP(MATVEC):', t_op, r_op/t_op, &
          & r_op
  Else
     Write (iot, FMT=600) 'time in OP:', t_op
  End If
  If (t_orth .Gt. zero) Then
     Write (iot, FMT=610) 'Re-Orthogonalization:', t_orth, &
          & r_orth/t_orth, r_orth
  Else
     Write (iot, FMT=600) 'time in orth:', t_orth
  End If
  If (t_res .Gt. zero) Then
     Write (iot, FMT=610) 'Restarting:', t_res, r_res/t_res,&
          & r_res
  Else
     Write (iot, FMT=600) 'time in restarting:', t_res
  End If
  If (t_tot .Gt. zero) Then
     Write (iot, FMT=610) 'TRLAN on this PE:', t_tot, r_tot/t_tot, &
          & r_tot
  Else
     Write (iot, FMT=600) 'total time in TRLAN:', t_tot
  End If
  If (info%verbose.Gt.0 .And. info%log_io.Ne.iot) Then
     Write(iot, *) 'Debug infomation written to files ',&
          & info%log_file(1:Index(info%log_file, ' ')-1), '####'
  Endif
  If (info%guess.Gt.1 .And. info%wrds_in.Gt.0) Then
     If (info%oldcpf .Ne. '') Then
        Write(iot,*) 'TRLAN restarted with checkpoint files ',&
             & info%oldcpf(1:Index(info%oldcpf,' ')-1), '####'
     Else
        Write(iot,*) 'TRLAN restarted with checkpoint files ',&
             & info%cpfile(1:Index(info%cpfile,' ')-1), '####'
     End If
     Write(iot,700) 'read', r_in, t_in, r_in/t_in
  End If
  If (info%clk_out.Gt.0 .And. info%wrds_out.Gt.0) Then
     Write(iot,*) 'Checkpoint files are ', &
          & info%cpfile(1:Index(info%cpfile, ' ')-1), '####'
     Write(iot,700) 'written', r_out, t_out, r_out/t_out
  End If
100 Format('Number of ', A, ' eigenpairs:', T32, I10, &
         &' (computed),', I11, ' (wanted)')
200 Format(A, T32, I10, ' (MAX:', I17, ' ).')
300 Format('Problem size: ', T32, I10, ' (PE:', I4, '),', I12,&
         & ' (Global)')
400 Format(A, T32, 1P, E10.3, ' (rel),', E16.3, ' (abs)')
500 Format(A, T32, I10)
600 Format(A, T25, 1P, E17.5, ' sec')
610 Format(A, T23, 1P, E12.4, ' sec,', E12.4, ' FLOP/S (', E11.4, ' FLOP)')
700 Format('Bytes ', A7, ': ', 1P, E12.5, ', Time(sec): ', E12.5,&
         & ', Rate(B/s): ', E12.5)
  If (info%npes .Gt. 1) Then
     ! write global performance information
     rinv = one / info%npes
     tmp1(1:12) = tmp1(1:12) * rinv
     Where (tmp1(1:6) .Gt. 0)
        tmp1(7:12) = tmp1(7:12) / tmp1(1:6)
     Elsewhere
        tmp1(7:12) = 0.0D0
     End Where
     If (tmp1(5).Eq.tmp1(6) .And. tmp1(5).Eq.zero) Then
        Write(iot, FMT=810)
        Write(iot, FMT=820) tmp1(1:4)
        Write(iot, FMT=830) tmp1(7:10)
     Else
        Write(iot, FMT=800)
        Write(iot, FMT=820) tmp1(1:6)
        Write(iot, FMT=830) tmp1(7:12)
     End If
  End If
800 Format(' -- Global summary --', /&
         &13X, 'Overall', 5X, 'MATVEC', 4X, 'Re-orth', 4X, 'Restart',&
         & 7X, 'Read', 6X, 'Write')
810 Format(' -- Global summary --', /&
         &13X, 'Overall', 5X, 'MATVEC', 4X, 'Re-orth', 4X, 'Restart')
820 Format('Time(ave)', 1P, 6E11.4)
830 Format('Rate(tot)', 1P, 6E11.4)
  If (iot.Eq.info%log_io .And. iot.Ne.6) Close(iot)
  Return
End Subroutine trl_print_info
!!!
! a more compact version of trl_print_info
!
! this is a local routine, indivadual PE can call it without regard of
! whether other PEs do the same and the output may be written to a
! different I/O unit number than the log_io
!!!
Subroutine trl_terse_info(info, iou)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T), Intent(in) :: info
  Integer, Intent(in) :: iou
  Integer :: rate
  Real(8) :: t_tot, t_op, t_orth, t_res
  If (info%clk_rate.Gt.0) Then
     t_op  = (info%tick_o+info%clk_op) / Dble(info%clk_rate)
     t_tot = (info%tick_t+info%clk_tot) / Dble(info%clk_rate)
     t_res  = (info%tick_r+info%clk_res) / Dble(info%clk_rate)
     t_orth = (info%tick_h+info%clk_orth) / Dble(info%clk_rate)
  Else
     ! get clock rate
     Call System_clock(count_rate=rate)
     t_op  = (info%tick_o+info%clk_op) / Dble(rate)
     t_tot = (info%tick_t+info%clk_tot) / Dble(rate)
     t_res  = (info%tick_r+info%clk_res) / Dble(rate)
     t_orth = (info%tick_h+info%clk_orth) / Dble(rate)
  Endif
  If (info%lohi.Gt.0) Then
     Write (iou, FMT=100) info%maxlan, info%restart, '+', info%ned,&
          & info%nec
  Else If (info%lohi.Lt.0) Then
     Write (iou, FMT=100) info%maxlan, info%restart, '-', info%ned,&
          & info%nec
  Else
     Write (iou, FMT=100) info%maxlan, info%restart, '0', info%ned,&
          & info%nec
  End If
  Write (iou, FMT=200) info%matvec, info%north, info%nloop, info%locked
  If (t_tot.Gt.1.0D-3 .And. Max(t_tot,t_op,t_res,t_orth).Lt.1.0D3) Then
     Write (iou, FMT=400) t_tot, t_op, t_orth, t_res
  Else
     Write (iou, FMT=300) t_tot, t_op, t_orth, t_res
  End If
100 Format('MAXLAN:', I10, ', Restart:', I10, ',   NED:', 2X, A1,&
         & I7, ',     NEC:', I10)
200 Format('MATVEC:', I10, ',  Reorth:', I10, ', Nloop:', I10,&
         & ', Nlocked:', I10)
300 Format('Ttotal:', 1PE10.3, ',    T_op:', 1PE10.3, ', Torth:',&
         & 1PE10.3, ',  Tstart:', 1PE10.3)
400 Format('Ttotal:', F10.6, ',    T_op:', F10.6, ', Torth:',&
         & F10.6, ',  Tstart:', F10.6)
  Return
End Subroutine trl_terse_info
!!!
! trl_check_ritz
! performs a standard check on the computed Ritz pairs
!
! Arguments:
! op    -- the operator, aka the matrix-vector multiplication routine
! nrow  -- the problem size
! rvec  -- the array storing the Ritz vectors
! alpha -- The ritz values
! beta  -- The residual norms returned from a Lanczos routine (optional)
! eval  -- The actual eigenvalues (optional)
! wrk   -- workspace (optional)
!!!
Subroutine trl_check_ritz(op, info, nrow, rvec, alpha, beta, eval, wrk)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T), Intent(in) :: info
  Integer, Intent(in) :: nrow
  Real(8), Dimension(:,:), Intent(in) :: rvec
  Real(8), Dimension(:), Intent(in) :: alpha
  Real(8), Dimension(:), Intent(in), Optional :: beta, eval
  Real(8), Dimension(:), Optional, Target :: wrk
  Interface
     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
     Subroutine op(nrow, ncol, xx, ldx, yy, ldy)
       Implicit None
       Integer, Intent(in) :: nrow, ncol, ldx, ldy
       Real(8), Intent(in) :: xx(ldx*ncol)
       Real(8), Intent(out) :: yy(ldy*ncol)
     End Subroutine op
     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
  End Interface
  ! local variables
  ! aq -- store the result of matrix-vector multiplication, size nrow
  ! rq -- store the Rayleigh-Quotient and the residual norms
  ! gsumwrk -- workspace left over for trl_g_sum to use
  Character(132) :: filename
  Integer :: i,lde,nev,lwrk,iaq,irq,iou
  Real(8) :: gapl, gapr
  Real(8), Dimension(:), Pointer :: aq, rq, gsumwrk, res, err
  ! dimension of the input arrays