trlaux.f90

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

! $Id: trlaux.f90,v 1.1 2004/05/04 06:29:35 fowlkes Exp $
!!!
! This file contains most auxillary routines which are not extensively
! used in during the normal operations of the TRLan, e.g., printing,
! error checking, etc.
!!!
! print an 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(1:), Intent(in) :: array
  If (Size(array).Gt.3) Then
     Write (info%log_io, *) 'PE', info%my_pe, ': ',&
          & Trim(title)
     Write (info%log_io, FMT=100) array
  Else
     Write (info%log_io, *) 'PE', info%my_pe, ': ',&
          & Trim(title), array
  End If
100 Format(8I10)
End Subroutine trl_print_int
!!!
! print a real(8) 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(1:), Intent(in) :: array
  If (Size(array).Gt.3) Then
     Write (info%log_io, *) 'PE', info%my_pe, ': ',&
          & Trim(title)
     Write (info%log_io, FMT=100) array
  Else If (Size(array).Gt.1) Then
     Write(info%log_io, FMT=*) 'PE', info%my_pe, ': ',&
          & Trim(title)
     Write (info%log_io, FMT=100) array
  Else
     Write(info%log_io, FMT=*) 'PE', info%my_pe, ': ',&
          & Trim(title), array
  End If
100 Format(1P,8E10.3)
End Subroutine trl_print_real
!!!
! current progress of eigenvalue solution
!!!
Subroutine trl_print_progress(info)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T), Intent(in) :: info
  Write(info%log_io, FMT=100) info%matvec, info%nloop, info%nec,&
       & info%north, info%nrand, info%stat
  Write(info%log_io, FMT=200) info%trgt, info%tres, info%crat
100 Format('MATVEC:', I10, ',    Nloop:', I10, ',     Nec:', I10, /,&
       & 'Reorth:', I10, ',    Nrand:', I10, ',    Ierr:', I10)
200 Format('Target:', 1P,E10.3, ',   ResNrm:', 1P,E10.3, ',   CFact:',&
         & 1P,E10.3)
End Subroutine trl_print_progress
!!!
! check orthogonality of the basis given
! -- used to debug of TRLANCZOS
!!!
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)
  Interface
     Subroutine trl_g_dot(mpicom, nrow, v1, ld1, m1, v2, ld2, m2, rr, wrk)
       Implicit None
       Integer, Intent(in) :: mpicom, nrow, ld1, ld2, m1, m2
       Real(8), Intent(in) :: v1(ld1,m1), v2(ld2,m2), rr(nrow)
       Real(8), Intent(out) :: wrk(m1+m2+m1+m2)
     End Subroutine trl_g_dot
  End Interface
  !
  ! local variables
  !
  Real(8), Parameter :: one=1.0D0, zero=0.0D0
  Integer :: i, j, jnd
  Real(8) :: nrmfro, nrminf
  !
  jnd = j1 + j2
  nrmfro = zero
  nrminf = zero
  If (jnd .Le. 0) Return
  If (lwrk .Lt. (jnd+jnd)) Then
     Write(info%log_io, *) 'TRL_CHECK_ORTH: workspace too small.'
     Return
  End If
  Write(info%log_io, *) 'TRL_CHECK_ORTH: check orthogonality of ',&
       & jnd, ' basis vectors.'
  !
  ! check orthognality of the basis vectors
  !
  Do i = 1, Min(j1, info%ntot) ! check no more than ntot vectors
     Call trl_g_dot(info%mpicom, nrow, v1, ld1, i, v2, ld2, 0, &
          & v1(1,i), wrk)
     wrk(i) = wrk(i) - one
     If (info%verbose.Gt.7) Then
        Write(info%log_io, *) 'Orthogonality level of v(', i,&
             & ') ..'
        Write(info%log_io, '(1P,8E10.3)') wrk(1:i)
     End If
     nrmfro = nrmfro + 2*Dot_product(wrk(1:i-1), wrk(1:i-1)) + &
          & wrk(i)*wrk(i)
     wrk(i+1) = Maxval(Abs(wrk(1:i)))
     nrminf = Max(nrminf, wrk(i+1))
  End Do
  Do i = 1, Min(j2, info%ntot-j1) ! check no more than ntot vectors
     j = j1 + i
     Call trl_g_dot(info%mpicom, nrow, v1, ld1, j1, v2, ld2, i,&
          & v2(1,i), wrk)
     wrk(j) = wrk(j) - one
     If (info%verbose.Gt.7) Then
        Write(info%log_io, *) 'Orthogonality level of v(', j, ') ..'
        Write(info%log_io, '(1P,8E10.3)') wrk(1:j)
     End If
     nrmfro = nrmfro + 2*Dot_product(wrk(1:j-1), wrk(1:j-1)) + &
          & wrk(j)*wrk(j)
     nrminf = Max(nrminf, Maxval(Abs(wrk(1:j))))
  End Do
  Write(info%log_io, FMT=100) info%matvec, jnd, Sqrt(nrmfro)
  Write(info%log_io, FMT=200) nrminf
100 Format('Frobenius norm of orthogonality level ', I10, I4, 1P, E14.5)
200 Format('Maximum abs. value of orthogonality level is ', 1P, E14.5)
End Subroutine trl_check_orth
!!!
! check Lanczos recurrence relation for debug purpose
!!!
Subroutine trl_check_recurrence(op, info, nrow, v1, ld1, m1, v2, ld2,&
     & m2, kept, alpha, beta, wrk, lwrk)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T), Intent(in) :: info
  Integer, Intent(in) :: nrow, ld1, ld2, m1, m2, kept, lwrk
  Real(8), Intent(in), Target :: v1(ld1,m1), v2(ld2,m2)
  Real(8), Intent(in) :: alpha(m1+m2-1), beta(m1+m2-1)
  Real(8), Target :: wrk(lwrk)
  External op
  Interface
     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(1:), Intent(in) :: array
     End Subroutine trl_print_real
  End Interface
  !
  ! local variables
  !
  Integer :: i, ii, j, j1, j2, jnd, mv1
  Character*80 :: title
  Real(8), Dimension(:), Pointer :: aq, qkp1, cs, alf, bet
  Real(8), Parameter :: zero=0.0D0, one=1.0d0
  !
  mv1 = m1
  If (m2 .Gt. 0) Then
     j2 = m2 -1
     j1 = m1
  Else
     j2 = 0
     j1 = m1 - 1
  End If
  jnd = j1 + j2
  If (lwrk .Lt. jnd*4+Max(jnd*4,nrow)) Then
     Write(info%log_io, *)&
          & 'TRL_CHECK_RECURRENCE: not enough workspace.'
     Return
  End If
  If (lwrk .Ge. jnd*4+nrow) Then
     aq => wrk(lwrk-nrow+1:lwrk)
  Else If (lwrk .Ge. jnd*4) Then
     Allocate(aq(nrow), stat=i)
     If (i.Ne.0) Then
        Write(info%log_io, *)&
             & 'TRL_CHECK_RECURRENCE: failed to allcoate workspace.'
        Return
     End If
  End If
  wrk(1:jnd*4) = zero
  cs => wrk(jnd+1:jnd+jnd)
  alf => wrk(jnd+jnd+1:jnd+jnd+jnd)
  bet => wrk(jnd+jnd+jnd+1:jnd*4)
  !
  ! first type of relation
  ! A q_i = Alpha_i q_i + Beta_i q_{k+1}
  !
  If (kept.Lt.Size(v1,2)) Then
     qkp1 => v1(1:nrow, kept+1)
  Else
     qkp1 => v2(1:nrow, kept-j1+1)
  End If
  Do i = 1, Min(j1, kept)
     Call op(nrow, 1, v1(1,i), nrow, aq, nrow)
     Do ii = 1, nrow
        alf(i) = alf(i) + aq(ii) * v1(ii,i)
        aq(ii) = aq(ii) - alpha(i) * v1(ii,i)
        bet(i) = bet(i) + aq(ii) * aq(ii)
        cs(i)  = cs(i)  + aq(ii) * qkp1(ii)
        aq(ii) = aq(ii) - beta(i) * qkp1(ii)
        wrk(i) = wrk(i) + aq(ii) * aq(ii)
     End Do
  End Do
  Do i = 1, kept-j1
     j = i + j1
     Call op(nrow, 1, v2(1,i), nrow, aq, nrow)
     Do ii = 1, nrow
        alf(j) = alf(j) + aq(ii) * v2(ii, i)
        aq(ii) = aq(ii) - alpha(j) * v2(ii, i)
        bet(j) = bet(j) + aq(ii) * aq(ii)
        cs(j)  = cs(j)  + aq(ii) * qkp1(ii)
        aq(ii) = aq(ii) - beta(j) * qkp1(ii)
        wrk(j) = wrk(j) + aq(ii) * aq(ii)
     End Do
  End Do
  !
  ! the (k+1)st base vector need to orthogonalize against all previous
  ! vectors
  !
  If (jnd .Gt. kept) Then
     Call op(nrow, 1, qkp1, nrow, aq, nrow)
     alf(kept+1) = Dot_product(aq, qkp1)
     aq = aq - alpha(kept+1) * qkp1
     Do i = 1, Min(j1,kept)
        aq = aq - beta(i) * v1(1:nrow,i)
     End Do
     Do i = 1, kept-j1
        j = j1 + i
        aq = aq - beta(j) * v2(1:nrow,i)
     End Do
     bet(kept+1) = Dot_product(aq, aq)
     If (kept+2 .Le. j1) Then
        cs(kept+1) = Dot_product(aq, v1(1:nrow, kept+2))
        aq = aq - beta(kept+1)*v1(1:nrow, kept+2)
     Else
        cs(kept+1) = Dot_product(aq, v2(1:nrow, kept+2-j1))
        aq = aq - beta(kept+1)*v2(1:nrow, kept+2-j1)
     End If
     wrk(kept+1) = Dot_product(aq, aq)
  End If
  !
  ! the third kind of relation -- normal three term recurrence
  ! depending the fact that if the lower-bound of loop is less than
  ! upper bound, the look should not be executed
  !
  Do i = kept+2, j1