trlan.f90

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

! $Id: trlan.f90,v 1.1 2004/05/04 06:29:35 fowlkes Exp $
!!! Top (user) level routines
!\Description:
! A thick-restart Lanczos routine for computing eigenvalues and
! eigenvectors of a real symmetric operator/matrix (A).
! -- only accept one input vector, the input vector is expected
!    to be stored in the (nec+1)st column of EVEC.
! -- it extends the Lanczos basis one vector at a time.
! -- orthogonality among the Lanczos vectors are maintained using
!    full re-orthogonalization when necessary.
!
!\Requirements:
! 1) User supplies OP with the specified interface
! 2) If (info%nec>0), evec(1:nrow, 1:info%nec) must contain valid
! eigenvectors and eval(1:nec) must be the corresponding eigenvalues.
! These eigenpairs are assumed to have zero residual norm inside TRLAN.
! 3) lde >= nrow
! 4) The arrays evec and eval must be large enough to store the
! solutions, i.e., mev >= info%ned and mev >= info%nec
! 5) The array wrk may be of arbitrary size.  Internally, the workspace
! size is
!
!        nrow*max(0,info%ned-size(evec,2))+maxlan*(maxlan+10)
!
! If wrk is smaller than this, trlan routine will allocate additional
! workspace to accommodate.
!
!\Arguments:
! op   -- the operator. Given a set of vectors X, it returns
!         op(X) == A*X
! info -- data structure to store the information about the eigenvalue
!         problem and the progress of TRLAN
! nrow -- number of rows that is on this processor (local problem size)
! mev  -- the number of Ritz pairs can be stored in eval and evec.
! eval -- array to store the eigenvalues
! evec -- array to store the eigenvectors
! lde  -- the leading dimension of the array evec (lde >= nrow)
! wrk  -- (optional) workspace, if it is provided and there is enough
!         space, the residual norm of the converged eigenpairs will be
!         stored at wrk(1:info%nec) on exit.
! lwrk -- (optional) the size of WRK.  When both WRK and LWRK are
!         present, than LWRK should correctly indicate the size of WRK.
!         If WRK is present by not LWRK, the size of WRK is assumed to
!         be MEV which is only enough to store the residual norms on
!         exit.  If WRK is not present, LWRK is not used even if it is
!         present.
!
! NOTE ON USING EXPLICIT SIZE ARRARY ARGUMENTS !
! This interface can be shorten to use assumed shape array arguments. 
! However, because it may cause extra copying of arrays when calling
! BLAS and LAPACK routines, we have decided to use only explicit size
! arrary arguments for most of the computing procedures.  Only the
! printing routines and debugging routines use assumed shape array
! arguments.
!
! the operator that defines the eigenvalue problem is expected to have
! the following interface
!  Subroutine OP(nrow, ncol, xin, ldx, yout, ldy)
!    Integer, Intent(in) :: nrow, ncol, ldx, ldy
!    Real(8), Dimension(ldx*ncol), Intent(in) :: xin
!    Real(8), Dimension(ldy*ncol), Intent(out) :: yout
!  End Subroutine OP
! Where the arguments are
! nrow -- number of rows in Xin and Yout on this PE (local problem
!         size)
! ncol -- number of columns in Xin and Yout, i.e., the number of
!         vectors to be multiplied
! xin  -- input vector to be multiplied.  It consists of Ncol column
!         vectors with each column stored in consecutive order
! ldx  -- the leading dimension of the array Xin.  The i-th column
!         vector starts with element (i-1)*ldx+1 and ends with element
!         (i-1)*ldx+nrow in Xin.
! yout -- the result array.  It stores the result of matrix-vector
!         multiplications.
! ldy  -- the leading dimension of the array yout.  The i-th column
!         vector starts with element (i-1)*ldy+1 in Yout and ends with
!         element (i-1)*ldy+nrow.
!
! Since it does not have an interface definition in Fortran 90, it is
! assumed to have an old Fortran 77 interface.  The user may declear
! the arrays Xin and Yout as two-dimension Fortran arrays in their
! implementation of the operator.
!!!
Subroutine trlan(op, info, nrow, mev, eval, evec, lde, wrk, lwrk)
  Use trl_info
  Implicit None
  External op
  Type(TRL_INFO_T) :: info
  Integer, Intent(in) :: nrow, mev, lde
  Real(8) :: eval(mev), evec(lde,mev)
  Integer, Intent(in), Optional :: lwrk
  Real(8), Target, Optional :: wrk(*)
  Interface
     !  ! the operator that defined the eigenvalue problem
     !  Subroutine OP(nrow, ncol, xin, ldx, yout, ldy)
     !    Integer, Intent(in) :: nrow, ncol, ldx, ldy
     !    Real(8), Dimension(ldx*ncol), Intent(in) :: xin
     !    Real(8), Dimension(ldy*ncol), Intent(out) :: yout
     !  End Subroutine OP
     ! the actual workhorse of the restarted Lanczos routine
     Subroutine trlanczos(op, info, nrow, mev, eval, evec, lde, base,&
          & ldb, nbas, wrk, lwrk)
       Use trl_info
       Implicit None
       Type(TRL_INFO_T) :: info
       Integer, Intent(in) :: nrow, mev, lde, ldb, nbas, lwrk
       Real(8) :: eval(1:mev)
       Real(8), Target :: evec(1:lde,1:mev), base(ldb*nbas), wrk(1:lwrk)
       External op
     End Subroutine trlanczos
     Subroutine trl_clear_counter(info, nrow, mev, lde)
       Use trl_info
       Implicit None 
       Type(TRL_INFO_T) :: info
       Integer, Intent(in) :: nrow, mev, lde
     End Subroutine trl_clear_counter
     Subroutine trl_print_setup(info, lbas, lmis, lwrk)
       Use trl_info
       Implicit None 
       Type(TRL_INFO_T), Intent(in) :: info
       Integer, Intent(in) :: lbas, lmis
       Integer, Intent(in), Optional :: lwrk
     End Subroutine trl_print_setup
     ! find the minimum flag value of 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
     Subroutine trl_time_stamp(iou)
       Implicit None
       Integer, Intent(in) :: iou
     End Subroutine trl_time_stamp
  End Interface
  Character(132) :: filename
  !
  ! The main task of this subroutine is prepare workspace required and
  ! check arguments
  !
  Integer :: ii, nbas, nmis, ibas, imis, ldb, lwrk0
  !
  ! local arrays -- all declared as 1-D array though they may
  ! accutually be used as 2-D arrays.
  !
  Real(8), Dimension(:), Pointer :: base, misc
  imis = -1    ! if this routine allocated misc, imis will be 0
  ibas = -1    ! if this routine allocated base, ibas will be 0
  Nullify(base, misc)
  Call System_clock(COUNT=ii, COUNT_RATE=info%clk_rate,&
       & COUNT_MAX=info%clk_max)
  info%clk_tot = ii
  !
  If (info%ned > mev) Then
     Print *, 'info%ned (', info%ned, ') is larger than mev (', mev,&
          & ')', 'reducing info%ned'
     info%ned = mev
  End If
  !
  ! there is nothing to do if there is no more eigenvalue to compute
  If (info%ned .Le. info%nec .Or. info%ned .Le. 0) Goto 888
  !
  ! determine the workspace size and leading dimensions
  If (Present(wrk) .And. Present(lwrk)) Then
     lwrk0 = lwrk
  Else If (Present(wrk)) Then
     lwrk0 = mev
  Else
     lwrk0 = 0
  End If
  info%stat = 0
  ldb = ((nrow+3)/4)*4
  If (Mod(ldb,4096) .Eq. 0) ldb = ldb + 4
  Call trl_clear_counter(info, nrow, mev, lde)
  If (info%stat .Ne. 0) Goto 888
  !
  ! Internally, the workspace is broken into two parts
  ! one to store (maxlan+1) Lanczos vectors, and the other to
  ! store all others (size maxlan*(maxlan+ned+14))
  ! The next If-block decides how the two arrays are mapped.
  !
  nbas = Max(1, info%maxlan-mev+1)
  ii = nbas * ldb
  nmis = info%maxlan*(info%maxlan+10)
  If (lwrk0 .Gt. Min(ii, nmis)) Then
     ! use wrk either as base or misc or both depending its size
     If (lwrk0 .Ge. ii+nmis) Then
        ! WRK is large enough for both arrays
        base => wrk(1:ii)
        misc => wrk(ii+1:lwrk0)
        nmis = lwrk0 - ii
     Else If (lwrk0 .Ge. Max(ii, nmis)) Then
        ! associate the larger one of base and misc to WRK
        If (ii .Ge. nmis) Then
           base => wrk(1:ii)
           Allocate(misc(nmis), stat=imis)
           If (imis .Ne. 0) info%stat = -4
        Else
           misc => wrk(1:lwrk0)
           nmis = lwrk0
           Allocate(base(ii), stat=ibas)
           If (ibas .Ne. 0) info%stat = -5
        End If
     Else If (ii .Le. nmis) Then
        ! base is smaller, associate base with WRK
        base => wrk(1:ii)
        Allocate(misc(nmis), stat=imis)
        If (imis .Ne. 0) info%stat = -4
     Else
        ! misc is smaller, associate misc with WRK
        misc => wrk(1:lwrk0)
        nmis = lwrk0
        Allocate(base(ii), stat=ibas)
        If (ibas .Ne. 0) info%stat = -5
     End If
  Else
     ! have to allocate both base and misc
     Allocate(base(ii), stat=ibas)
     If (ibas .Ne. 0) info%stat = -5
     Allocate(misc(nmis), stat=imis)
     If (imis .Ne. 0) info%stat = -4
  End If
  !
  ! make sure every process is successful so far
  ii = trl_sync_flag(info%mpicom, info%stat)
  info%stat = ii
  If (ii .Ne. 0) Goto 888
  !
  ! open debug files if desired
  If (info%verbose.Gt.0 .And. info%log_io.Ne.6) Then
     filename = trl_pe_filename(info%log_file, info%my_pe, info%npes)
     Open(info%log_io, file=filename, status='replace',&
          & action='write', position='rewind', iostat=ii)
     If (ii .Ne. 0) info%log_io = 6
  End If
  If (info%verbose .Gt. 0) Then
     Call trl_time_stamp(info%log_io)
     Call trl_print_setup(info, nbas*ldb, nmis, lwrk0)
  End If
  !
  ! call trlanczos to do the real work
  base = 0.0D0
  misc = 0.0D0
  Call trlanczos(op, info, nrow, mev, eval, evec, lde, base, ldb, nbas,&
       & misc, nmis)
  If (info%verbose.Gt.0 .And. info%log_io.Ne.6) Close(info%log_io)
  ii = Max(info%nec, info%ned)
  If (lwrk0 .Ge. ii) Then
     wrk(1:ii) = misc(1:ii)
  End If
  !
  ! DONE, reclaim the space allocated
888 If (imis .Eq. 0) Deallocate(misc)
  If (ibas .Eq. 0) Deallocate(base)
  Call System_clock(count=ii)
  If (ii .Ge. info%clk_tot) Then
     info%clk_tot = ii - info%clk_tot
  Else
     info%clk_tot = ii + (info%clk_max - info%clk_tot)
  End If
  Return
End Subroutine trlan
!!!
! clear the counters inside info and performs a minimal check on the
! input parameters
!!!
Subroutine trl_clear_counter(info, nrow, mev, lde)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T) :: info
  Integer, Intent(in) :: nrow, mev, lde
  Integer, External :: trl_sync_flag
  Integer :: ntmp
  !!
  info%stat = 0
  If (nrow.Ne.info%nloc .Or. nrow.Gt.info%ntot) Then
     Print *, 'TRLAN: ''info'' not setup for this problem.'
     Print *, '       Please reset ''info'' before calling TRLAN.'
     info%stat = -1
  End If
  If (info%ned+info%ned >= info%ntot) Then
     Print *, 'TRLAN: info%ned (', info%ned, ') is large relative ', &
          & 'to the matrix dimension (', info%ntot, ')'
     Print *, '       it might be more appropriate to use LAPACK ', &
          & 'dsyev/ssyev'
  End If
  If (info%nec .Lt. 0) info%nec = 0
  If (lde .Lt. nrow) Then
     Print *, 'TRLAN: leading dimension of EVEC to small.'
     info%stat = -2
  End If
  If (.Not. (info%tol .Lt. 1D0)) Then
     info%tol = Sqrt(Epsilon(info%tol))
  Else If (.Not. (info%tol .Gt. Tiny(info%tol))) Then
     info%tol = Epsilon(info%tol)
  End If
  info%ned = Min(info%ntot, info%ned)
  If (mev .Lt. info%ned) Then
     Print *, 'TRLAN: array EVAL and EVEC can not hold wanted', &
          & ' number of eigenpairs.'
     info%stat = -3
  End If
  info%nrand = info%stat
  info%stat = trl_sync_flag(info%mpicom, info%nrand)
  !
  ! decide what is a good maximum basis size to use
  !
  If (info%maxlan .Le. Min(info%ned+3, info%ntot)) Then
     info%maxlan = info%ned+Min(info%ned,20)+Int(Log(Dble(info%ntot)))
     info%maxlan = Min(info%maxlan, info%ntot)
  End If
  If (info%maxlan < mev) Then
     ntmp = Min(info%ntot, Max(100+info%ned, 10*info%ned))
     If (mev < ntmp) Then
        info%maxlan = mev
     Else
        info%maxlan = ntmp;
     End If
  End If
  !!! clear regular counters
  info%tmv = -1
  info%klan = Min(info%maxlan, info%ntot)
  info%locked = info%nec
  info%matvec = 0
  info%nloop = 0
  info%north = 0
  info%nrand = 0
  info%flop = 0
  info%rflp = 0
  info%flop_h = 0
  info%rflp_h = 0
  info%flop_r = 0
  info%rflp_r = 0
  info%tick_t = 0.0D0
  info%clk_op = 0
  info%tick_o = 0.0D0
  info%clk_orth = 0
  info%tick_h = 0.0D0
  info%clk_res = 0
  info%tick_r = 0.0D0
  info%clk_in = 0
  info%clk_out = 0
  info%wrds_in = 0
  info%wrds_out = 0
  Return
End Subroutine trl_clear_counter
!!!
! print the definition of the eigenvalue proble
!!!
Subroutine trl_print_setup(info, lbas, lmis, lwrk)
  Use trl_info
  Implicit None 
  Type(TRL_INFO_T), Intent(in) :: info
  Integer, Intent(in) :: lbas, lmis
  Integer, Intent(in), Optional :: lwrk
  !
  ! print the problem parameters
  If (info%lohi .Gt. 0) Then
     Write (info%log_io, FMT=100) info%ned, 'largest'
  Else If (info%lohi .Lt. 0) Then
     Write (info%log_io, FMT=100) info%ned, 'smallest'
  Else
     Write (info%log_io, FMT=100) info%ned, 'first converged'
  Endif
  Write(info%log_io, FMT=150) info%nloc, info%my_pe, info%ntot
  Write(info%log_io, FMT=200) 'Maximum basis size:', info%maxlan
  Write(info%log_io, FMT=200) 'Dynamic restarting scheme:', info%restart
  Write(info%log_io, FMT=200) 'Maximum applications of the operator:',&
       & info%maxmv
  Write(info%log_io, FMT=300) 'Relative convergence tolerance:', info%tol
100 Format('TRLAN is to compute ', I6, 1X, A, ' eigenpair(s).')
150 Format('Problem dimension: ', I9, '(PE:', I4, '),', I12, '(Global)')
200 Format(A, T40, I10)
300 Format(A, T40, 1PE10.3)
  ! initial guess
  If (info%guess .Eq. 1) Then
     Write(info%log_io, *) 'User provided the starting vector.'
  Else If (info%guess .Eq. 0) Then
     Write(info%log_io, *) 'TRLAN uses [1,1,...] as starting vctor.'
  Else If (info%guess .Lt. 0) Then
     Write(info%log_io, *) 'TRLAN generates a random starting vector.'
  Else If (info%oldcpf .Ne. '') Then
     Write(info%log_io, *)&
          & 'Restarting with existing checkpoint files ',&
          & Trim(info%oldcpf), '####'
  Else
     Write(info%log_io, *)&
          & 'Restarting with existing checkpoint files ',&
          & Trim(info%cpfile), '####'
  End If
  If (info%cpflag .Gt. 0) Then
     Write(info%log_io, *) 'TLRAN will write about ', info%cpflag,&
          &' sets of checkpointing files ', Trim(info%cpfile), '####.'
  End If
  !
  ! print the workspace size parameters
  Write(info%log_io, *) '(required) array BASE size is ', lbas
  Write(info%log_io, *) '(required) array MISC size is ', lmis
  If (Present(lwrk)) Then
     If (lwrk .Gt. 0) Then
        Write(info%log_io, *) 'Caller has supplied a work array with ',&
          & lwrk, ' elements.'
     Else
        Write(info%log_io, *) 'Caller did not supply work array.'
     End If
  Else
     Write(info%log_io, *) 'Caller did not supply work array.'
  End If
End Subroutine trl_print_setup
!!!
! set information related to debugging, the initialization routine
! trl_init_info sets the parameters to not allow no debug information
! to be printed.
!
! arguments:
! info   -- the info to be modified
! msglvl -- the new message level
!           < 0 : nothing printed
!           1:10 -- the high the level, the more debug information is
!                printed
! iou    -- Fortran I/O unit to be used to write out the debug
!           information
! file   -- leading part of the debug file name.