trlcore.f90

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

! 1. if (global re-orthogonalization is needed)
!      call trl_cgs
!    else
!      perform extended local re-reorthogonalization
!    endif
! 2. perform normalization
!
! workspace requirement: lwrk >= 2*(m1+m2)
!!!
Subroutine trl_orth(nrow, v1, ld1, m1, v2, ld2, m2, rr, kept, alpha,&
     & beta, wrk, lwrk, info)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T) :: info
  Integer, Intent(in) :: nrow, ld1, ld2, m1, m2, kept, lwrk
  Real(8), Intent(in), Target :: v1(ld1,m1), v2(ld2,m2)
  Real(8) :: rr(nrow), alpha(m1+m2), beta(m1+m2), wrk(lwrk)
  !
  ! local variables
  !
  Real(8), Parameter :: zero=0.0D0, one=1.0D0
  Integer :: i, jnd, jm1, no, nr
  Real(8) :: tmp
  Real(8), Dimension(:), Pointer :: qa, qb
  Interface
     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
     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
  !
  ! check for workspace size
  !
  jnd = m1 + m2
  jm1 = jnd - 1
  tmp = zero
  If (ld1.Ge.nrow .And. ld2.Ge.nrow .And. lwrk.Ge.Max(4,jnd+jnd)) Then
     info%stat = 0
  Else
     info%stat = -101
     Return
  End If
!!$Print *, 'Entering TRL_ORTH', m1, m2 
!!$Do i = 1, m1
!!$Print *, 'v(:,',i, ') '
!!$Print *, v1(1:nrow,i)
!!$End Do
!!$Do i = 1, m2
!!$Print *, 'v(:,',m1+i,') '
!!$Print *, v2(1:nrow, i)
!!$End Do
!!$Print *, 'rr '
!!$Print *, rr
  !
  ! compute the norm of the vector RR
  !
  wrk(1) = Dot_product(rr, rr)
  Call trl_g_sum(info%mpicom, 1, wrk(1), wrk(2))
  If (.Not.(wrk(1).Ge.zero) .Or. .Not.(wrk(1).Le.Huge(tmp))) Then
     info%stat = -102
     Return
  End If
  beta(jnd) = Sqrt(wrk(1))
  tmp = alpha(jnd)*alpha(jnd)
  If (jm1 .Gt. kept) Then
     tmp = tmp + beta(jm1)*beta(jm1)
     info%flop = info%flop + 2*nrow + 4
  Else If (kept.Gt.0) Then
     tmp = tmp + Dot_product(beta(1:jm1), beta(1:jm1))
     info%flop = info%flop + 2*(nrow + kept + 2)
  End If
  !
  ! decide whether to perform global or local reothogonalization
  If (jnd .Ge. info%ntot) Then
     ! very specfial case: do nothing more
     info%stat = 0
     beta(jnd) = zero
  Else If (tmp.Ge.wrk(1) .Or. wrk(1).Le.Tiny(tmp) .Or. &
       & beta(jm1).Le.Tiny(tmp) .Or. jm1.Eq.kept) Then
     ! perform global re-orthogonalization
     nr = info%nrand
     no = info%north
     Call trl_cgs(info%mpicom, info%my_pe, nrow, v1, ld1, m1, v2, &
          & ld2, m2, rr, beta(jnd), alpha(jnd), info%north, &
          & info%nrand, wrk, info%stat)
     info%flop = info%flop + 4*nrow*((info%north-no)*jnd +&
          & info%nrand-nr) + nrow
  Else If (jnd .Gt. 1) Then
     ! perform local re-orthogonalization against two previous vectors
     If (m2.Gt.1) Then
        qa => v2(1:nrow, m2)
        qb => v2(1:nrow, m2-1)
     Else If (m2.Eq.1) Then
        qa => v2(1:nrow, 1)
        qb => v1(1:nrow, m1)
     Else
        qa => v1(1:nrow, m1)
        qb => v1(1:nrow, jm1)
     End If
     wrk(1) = zero
     wrk(2) = zero
     Do i = 1, nrow
        wrk(1) = wrk(1) + qa(i)*rr(i)
        wrk(2) = wrk(2) + qb(i)*rr(i)
     End Do
     Call trl_g_sum(info%mpicom, 2, wrk(1:2), wrk(3:4))
     alpha(jnd) = alpha(jnd) + wrk(1)
     rr = rr - wrk(1)*qa - wrk(2)*qb
     tmp = one / beta(jnd)
     rr = tmp*rr
     info%flop = info%flop + 9*nrow
  Else
     ! perform local re-orthogonalization against the only vector
     If (m1.Eq.1) Then
        qa => v1(1:nrow, 1)
     Else
        qa => v2(1:nrow, 1)
     End If
     wrk(1) = Dot_product(qa, rr)
     Call trl_g_sum(info%mpicom, 1, wrk(1:1), wrk(2:2))
     alpha(jnd) = alpha(jnd) + wrk(1)
     rr = rr - wrk(1)*qa
     tmp = one / beta(jnd)
     rr = tmp*rr
     info%flop = info%flop + 5*nrow
  End If
  ! when beta(jnd) is exceedingly small, it should be treated as zero
  If (info%stat .Eq. 0) Then
     If (beta(jnd) .Le. Epsilon(tmp)*Abs(alpha(jnd))) Then
        beta(jnd) = zero
     Else If (jnd .Ge. info%ntot) Then
        beta(jnd) = zero
     End If
  End If
End Subroutine trl_orth
!!!
! transforms an real symmetric arrow matrix into a
! symmetric tridiagonal matrix
!!!
Subroutine trl_tridiag(nd, alpha, beta, rot, alfrot, betrot, wrk, lwrk,&
     & ierr)
  Implicit None
  Integer, Intent(in) :: nd, lwrk
  Integer, Intent(out) :: ierr
  Real(8), Intent(in), Dimension(nd) :: alpha, beta
  Real(8) :: rot(nd*nd), alfrot(nd), betrot(nd), wrk(lwrk)
  External dsytrd, dorgtr
  !  Interface
  !     Subroutine dsytrd(uplo, n, a, lda, d, e, tau, wrk, lwrk, ierr)
  !       Character, Intent(in) :: uplo
  !       Integer, Intent(in) :: n, lda, lwrk
  !       Integer :: ierr
  !       Real(8) :: a(lda,n), d(n), e(n), tau(n), wrk(lwrk)
  !     End Subroutine dsytrd
  !     Subroutine dorgtr(uplo, n, a, lda, tau, wrk, lwrk, ierr)
  !       Character, Intent(in) :: uplo
  !       Integer, Intent(in) :: n, lda, lwrk
  !       Integer :: ierr
  !       Real(8) :: a(lda,n), tau(n), wrk(lwrk)
  !     End Subroutine dorgtr
  !  End Interface
  !
  ! local variables
  !
  Character, Parameter :: upper='U'
  Integer :: lwrk2
  !
  ! special case
  !
  If (nd .Le. 1) Then
     rot=1.0D0
     alfrot = alpha
     betrot = beta
     ierr = 0
     Return
  End If
  If (lwrk.Ge.nd+nd) Then
     ierr = 0
  Else
     ierr = -111
     Return
  End If
  !
  ! first form the array matrix as a full matrix in rot
  !
  rot = 0.0D0
  rot(1:nd*nd:nd+1) = alpha(1:nd)
  rot((nd-1)*nd+1:nd*nd-1) = beta(1:nd-1)
  rot(nd:nd*(nd-1):nd) = beta(1:nd-1)
  lwrk2 = lwrk - nd
  !
  ! call LAPACK routines to reduce the matrix into tridiagonal form
  ! and generate the rotation matrix
  !
  Call dsytrd(upper, nd, rot, nd, alfrot, betrot, wrk, wrk(nd+1), &
       &  lwrk2, ierr)
  If (ierr .Ne. 0) Then
     ierr = -112
     Return
  End If
  betrot(nd) = beta(nd)
  Call dorgtr(upper, nd, rot, nd, wrk, wrk(nd+1), lwrk2, ierr)
  If (ierr .Ne. 0) Then
     ierr = -113
     Return
  End If
End Subroutine trl_tridiag
!!!
! evaluates the eigenvalues and their corresponding residual norms of a
! real symmetric tridiagonal matrix
!
! it returns eigenvalues in two sections
! the first section is the locked eigenvalues, their residual norms are
! zero
! the second section contains the new Ritz values, they are in
! ascending order.
! res will contain corresponding residual norms
!!!
Subroutine trl_get_eval(nd, locked, alpha, beta, lambda, res, wrk,&
     & lwrk, ierr)
  Implicit None
  Integer, Intent(in) :: nd, locked, lwrk
  Integer :: ierr
  Real(8), Intent(in) :: alpha(nd), beta(nd)
  Real(8), Intent(out) :: lambda(nd), res(nd), wrk(lwrk)
  External dstqrb
  !  Interface
  !     Subroutine dstqrb(n, d, e, z, wrk, ierr)
  !       Integer, Intent(in) :: n
  !       Integer :: ierr
  !       Real(8) :: d(n), e(n), z(n), wrk(n+n-2)
  !     End Subroutine dstqrb
  !  End Interface
  If (lwrk .Ge. 3*nd) Then
     ierr = 0
  Else
     ierr = -121
     Return
  End If
  lambda = alpha
  wrk(1:nd-locked) = beta(locked+1:nd)
  Call dstqrb(nd-locked, lambda(locked+1:nd), wrk(1:nd-locked),&
       & res(locked+1:nd), wrk(1+nd:), ierr)
  If (ierr .Eq. 0) Then
     res(1:locked) = 0.0D0
     res(locked+1:nd) = beta(nd) * Abs(res(locked+1:nd))
  Else
     ierr = -122
  End If
End Subroutine trl_get_eval
!!!
! count the numer of wanted eigenvalues that have small residual norms
! eigenvalues are assumed to be order from small to large
!!!
Subroutine trl_convergence_test(nd, lambda, res, info, wrk)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T) :: info
  Integer, Intent(in) :: nd
  Real(8), Intent(in) :: lambda(nd), res(nd)
  Real(8), Intent(out) :: wrk(nd+nd)
  Real(8) :: bnd
  Integer :: i, j, ncl, ncr
  External dsort2
  !
  ! copy lambda and res to wrk, sort them in ascending order of lambda
  wrk(nd+1:nd+nd) = lambda(1:nd)
  wrk(1:nd) = Abs(res(1:nd))
  Call dsort2(nd, wrk(nd+1:nd+nd), wrk(1:nd))
  !
  ! determine the convergence rate of the previous target
  If (info%tmv.Gt.0 .And. info%matvec.Gt.info%tmv) Then
     j = 1
     bnd = Abs(lambda(j)-info%trgt)
     Do i = 1, nd
        If (Abs(lambda(i)-info%trgt) .Lt. bnd) Then
           bnd = Abs(lambda(i)-info%trgt)
           j = i
        End If
     End Do
     If (info%tres .Gt. res(j)) Then
        bnd = res(j) / info%tres
        If (bnd .Gt. 0.0D0) Then
           info%crat = Log(bnd)/Dble(info%matvec-info%tmv)
        Else
           info%crat = 0.0D0
        End If
     Else
        info%crat = 1.0D0
     End If
  End If
  !
  ! find out who has converged at the lower end of the spectrum
  info%anrm = Max(info%anrm, Abs(wrk(nd+1)), Abs(wrk(nd+nd)))
  bnd = Tiny(info%anrm) + info%tol*info%anrm
  ncl = 0
  ncr = nd+1
  If (info%lohi .Le. 0) Then
     ncl = nd
     i = 1
     Do While (i.Le.nd)
        If (wrk(i) .Lt. bnd) Then
           i = i + 1
        Else
           ncl = i - 1
           i = nd + 1
        End If
     End Do
  End If
  ! find out who has converged at the high end of the spectrum
  If (info%lohi .Ge. 0) Then
     ncr = 1
     i = nd
     Do While (i.Ge.1)
        If (wrk(i) .Lt. bnd) Then
           i = i - 1
        Else
           ncr = i + 1
           i = 0
        End If
     End Do
  End If
  ! determine the number of wanted eigenvalues that have converged
  ! compute the next target
  info%tmv = info%matvec
  If (info%lohi .Lt. 0) Then
     info%nec = ncl
     info%trgt = wrk(nd+Min(nd,ncl+1))
     info%tres = wrk(Min(nd,ncl+1))
  Else If (info%lohi .Gt. 0) Then
     info%nec = nd-ncr+1
     info%trgt = wrk(nd+Max(1,ncr-1))
     info%tres = wrk(Max(1,ncr-1))
  Else
     If (ncr .Le. ncl) Then
        ncl = nd/2
        ncr = ncl+1
        info%trgt = wrk(nd+(nd+1)/2)
        info%tres = wrk((nd+1)/2)
     Else If (wrk(ncl+1).Le.wrk(ncr-1)) Then
        info%trgt = wrk(nd+ncl+1)
        info%tres = wrk(ncl+1)
     Else
        info%trgt = wrk(nd+ncr-1)
        info%tres = wrk(ncr-1)
     End If
     info%nec = ncl + nd - ncr + 1
     Do i = ncl+1, ncr-1
        If (wrk(i) .Lt. bnd) info%nec = info%nec + 1
     End Do
  End If
End Subroutine trl_convergence_test
!!!
! sort the eigenvalues so that the wanted eigenvalues are ouputed to the
! user in the front of the arrays
! the final Ritz values are in ascending order so that
! DSTEIN can be used to compute the eigenvectors
!!!
Subroutine trl_sort_eig(nd, lohi, nec, lambda, res)
  Implicit None
  Integer, Intent(in) :: nd, lohi, nec
  Real(8) :: lambda(nd), res(nd)
  External dsort2, dsort2a
  Integer :: i, j
  If (lohi .Eq. 0) Then
     ! sort the eigenvalues according to their absolute values
     Call dsort2a(nd, lambda, res)
     ! sort the first nec eigenvalue in the order of lambda
     Call dsort2(nec, lambda, res)
  Else
     ! sort the eigenvalues and residual norms in ascending order of the
     ! eigenvalues
     Call dsort2(nd, lambda, res)
     If (lohi .Gt. 0) Then
        ! move the largest ones to the front (still ascending order)
        j = nd - nec + 1
        Do i = 1, nec