trlan.f90

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

  lde = Size(rvec, 1)
  nev = Size(rvec, 2)
  If (nev .Le. 0) Return  ! there is nothing to do
  If (nev .Gt. Size(alpha,1)) Then
     Print *, 'TRL_CHECK_RITZ: Rvec and Alpha array size do not match.'
     Return
  End If
  ! figure out whether it is necessary to allocate new workspace
  iaq = nrow
  irq = nev+nev+nev
  If (Present(wrk)) Then
     ! WRK provide by caller, try to use it
     lwrk = Size(wrk)
     If (lwrk .Ge. iaq+irq) Then
        aq => wrk(1:iaq)
        rq => wrk(iaq+1:iaq+irq)
        gsumwrk => wrk(iaq+irq+1:lwrk)
     Else If (lwrk .Ge. iaq) Then
        aq => wrk(1:iaq)
        gsumwrk => wrk(iaq+1:lwrk)
        Allocate(rq(nev+nev+nev), stat=irq)
        If (irq .Ne. 0) Then
           Print *, 'TRL_CHECK_RITZ: Failed to allocated workspace RQ.'
           Return
        End If
     Else If (lwrk .Ge. irq) Then
        rq => wrk(1:irq)
        gsumwrk => wrk(irq+1:lwrk)
        Allocate(aq(nrow), stat=iaq)
        If (iaq .Ne. 0) Then
           Print *, 'TRL_CHECK_RITZ: Failed to allocated workspace AQ.'
           Return
        End If
     Else
        gsumwrk => wrk
        Allocate(aq(nrow),  stat=iaq)
        If (iaq .Ne. 0) Then
           Print *, 'TRL_CHECK_RITZ: Failed to allocated workspace AQ.'
           Return
        End If
        Allocate(rq(nev+nev+nev), stat=irq)
        If (irq .Ne. 0) Then
           Print *, 'TRL_CHECK_RITZ: Failed to allocated workspace RQ.'
           Deallocate(aq)
           Return
        End If
     End If
  Else
     ! WRK not provided -- allocate space for AQ and RQ,
     ! gsumwrk points to the last third of RQ
     Allocate(aq(nrow),  stat=iaq)
     If (iaq .Ne. 0) Then
        Print *, 'TRL_CHECK_RITZ: Failed to allocated workspace AQ.'
        Return
     End If
     Allocate(rq(nev+nev+nev), stat=irq)
     If (irq .Ne. 0) Then
        Print *, 'TRL_CHECK_RITZ: Failed to allocated workspace RQ.'
        Deallocate(aq)
        Return
     End If
     gsumwrk => rq(nev+nev+1:nev+nev+nev)
  End If
  aq = 0
  rq = 0
  gsumwrk = 0
  !
  ! go through each Ritz pair one at a time, compute Rayleigh
  ! quotient and the corresponding residual norm
  res => rq(nev+1:nev+nev)
  Do i = 1, nev
     Call op(nrow, 1, rvec(:,i), lde, aq, nrow)
     ! Rayleigh quotient -- assuming rvec(:,i) has unit norm
     rq(i) = Dot_product(rvec(1:nrow,i), aq)
     Call trl_g_sum(info%mpicom, 1, rq(i:i), gsumwrk)
     aq = aq - rq(i)*rvec(1:nrow,i)
     res(i) = Dot_product(aq, aq)
  End Do
  Call trl_g_sum(info%mpicom, nev, rq(nev+1:nev+nev), gsumwrk)
  !
  ! compute the error estimate based on computed residual norms and the
  ! Ritz values
  err => rq(nev+nev+1:nev+nev+nev)
  res = Sqrt(res)
  gapl = alpha(nev) - alpha(1)
  Do i = 1, nev-1
     gapr = alpha(i+1) - alpha(i)
     gapl = Min(gapl, gapr)
     If (res(i).Ge.gapl) Then
        err(i) = res(i)
     Else
        err(i) = res(i) * res(i) / gapl
     End If
     gapl = gapr
  End Do
  If (res(nev).Ge.gapl) Then
     err(nev) = res(nev)
  Else
     err(nev) = res(nev) * res(nev) / gapl
  End If
  !
  ! if a log file was opened before, open it again to append the new
  ! stuff to the end
  iou = info%log_io
  If (iou.Le.0 .Or. info%verbose.Le.0) iou = 6
  If (iou .Ne. 6) Then
     filename = trl_pe_filename(info%log_file, info%my_pe, info%npes)
     Open(iou, file=filename, status='OLD', position='APPEND',&
          & action='WRITE', iostat=i)
     If (i .Ne. 0) iou = 6
  End If
  ! if writing to IO unit 6, only PE 0 does it
  If (iou.Eq.6 .And. info%my_pe.Gt.0) go to 888
  ! print out the information
  Write (iou, FMT=100)
  If (Present(beta) .And. Present(eval)) Then
     Do i = 1, nev
        Write(iou, FMT=200) alpha(i), res(i), beta(i)-res(i),&
             & err(i), rq(i)-alpha(i), eval(i)-alpha(i)
     End Do
  Else If (Present(beta)) Then
     Do i = 1, nev
        Write(iou, FMT=200) alpha(i), res(i),&
             & beta(i)-res(i), err(i), rq(i)-alpha(i)
     End Do
  Else If (Present(eval)) Then
     Do i = 1, nev
        Write(iou, FMT=300) alpha(i), res(i), err(i), rq(i)-alpha(i),&
             & eval(i)-alpha(i)
     End Do
  Else
     Do i = 1, nev
        Write(iou, FMT=300) alpha(i), res(i), err(i), rq(i)-alpha(i)
     End Do
  End If
100 Format('TRL_CHECK_RITZ:', /,&
         & 11X, 'Ritz value', 7X, 'res norm', '   res diff',&
         & '  est error', '  diff w rq', ' act. error')
200 Format(1P, G25.17, 5E11.3)
300 Format(1P, G25.17, E11.3, 11X, 3E11.3)
  If (iou.Ne.6) Close(iou)
  ! deallocate workspace allocated in this routine
888 If (irq .Eq. 0) Deallocate(rq)
  If (iaq .Eq. 0) Deallocate(aq)
End Subroutine trl_check_ritz
!!!
! A separate Rayleigh-Ritz projection routine
! Given a set of approximately orthonormal vectors (V), this routine
! performs the following operations
! (1) V'*V ==> G
! (2) R'*R :=  G
! (3) V'*A*V => H1, inv(R')*H1*inv(R) => H
! (4) Y*D*Y' := H
! (5) V*inv(R)*Y => V, diag(D) => lambda,
!     r(i) = ||A*V(:,i)-lambda(i)*V(:,i)||
!
! Arguments:
! op   -- the operator (matrix-vector multiplication routine)
! info -- the structure for storing information about TRLAN
! evec -- the Ritz vectors to be manipulated
! eres -- the array to store new Ritz values and residual norms
! base -- the (optional) workspace to store the result of MATVEC
! wrk  -- the (optional) workspace to store projection matrix, etc..
!!!
Subroutine trl_ritz_projection(op, info, evec, eres, wrk, base)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T) :: info
  Real(8), Dimension(:), Optional, Target :: wrk, base
  Real(8), Dimension(:, :), Target :: evec
  Real(8), Dimension(:) :: eres
  !! local variables
  Real(8), Parameter :: one=1.0D0, zero=0.0D0
  Integer :: i, ierr, lde, lwrk, mev, nev, nsqr, nrow, iuau, irvv
  Real(8), Dimension(:), Pointer :: rvv, uau, wrk2, avec
  External dgemm, dpotrf, dstev, dtrtrs
  Interface
     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
     Subroutine trl_ritz_vectors(nrow, lck, ny, yy, ldy, vec1, ld1, m1,&
          & vec2, ld2, m2, wrk, lwrk)
       Implicit None
       Integer, Intent(in) :: nrow, ny, lck, ldy, ld1, m1, m2, ld2, lwrk
       Real(8), Intent(in) :: yy(ldy, ny)
       Real(8) :: vec1(ld1,m1), vec2(ld2,m2), wrk(lwrk)
     End Subroutine trl_ritz_vectors
  End Interface
  !
  lde = Size(evec,1)
  mev = Size(evec,2)
  nrow = info%nloc
  If (info%nec .Gt. 0) Then
     nev = info%nec + 1
  Else
     nev = Min(info%ned, mev-1)
     If (info%lohi.Ne.0) nev = nev + 1
  End If
  nsqr = nev*nev
  If (Present(wrk)) Then
     lwrk = Size(wrk, 1)
  Else
     lwrk = 0
  End If
  If (Present(base)) Then
     avec => base(1:nrow)
  Else If (mev .Gt. nev) Then
     avec => evec(1:nrow, mev)
  Else
     Allocate(avec(nrow))
  End If
  ! memory allocation -- need 3*nev*nev elements, will allocate them
  ! in two consecutive blocks, uau(nev*nev), rvv(2*nev*nev)
  ! in actual use, rvv is further split in two until the last operation
  iuau = nsqr
  irvv = nsqr+nsqr
  If (lwrk .Ge. iuau+irvv) Then
     uau => wrk(1:nsqr)
     rvv => wrk(nsqr+1:nsqr+nsqr)
     wrk2 => wrk(nsqr+nsqr+1:lwrk)
  Else If (lwrk .Ge. irvv) Then
     rvv => wrk(1:nsqr)
     wrk2 => wrk(nsqr+1:lwrk)
     Allocate(uau(nsqr), stat=iuau)
     If (iuau .Ne. 0) Then
        info%stat = -231
        Goto 888
     End If
  Else If (lwrk .Ge. iuau) Then
     uau => wrk(1:nsqr)
     Allocate(rvv(nsqr+nsqr), stat=irvv)
     If (irvv .Ne. 0) Then
        info%stat = -232
        Goto 888
     End If
     wrk2 => rvv(nsqr+1:nsqr+nsqr)
  Else
     Allocate(uau(nsqr), stat=iuau)
     If (iuau .Ne. 0) Then
        info%stat = -231
        Goto 888
     End If
     Allocate(rvv(nsqr+nsqr), stat=irvv)
     If (irvv .Ne. 0) Then
        info%stat = -232
        Goto 888
     End If
     wrk2 => rvv(nsqr+1:nsqr+nsqr)
  End If
  ! step (1) : V'*V ==> G
  Call dgemm('T', 'N', nev, nev, nrow, one, evec, lde, evec, lde, zero,&
       & rvv, nev)
  Call trl_g_sum(info%mpicom, nsqr, rvv(1:nsqr), wrk2)
  ! step (2) : Choleskey factorization of G
  Call dpotrf('U', nev, rvv, nev, ierr)
  If (ierr .Ne. 0) Then
     info%stat = -234
     Goto 888
  End If
  ! step (3) : compute H_1 = V'*A*V
  ! use the first nrow elements of avec to store the results of
  ! matrix-vector multiplication
  wrk2 = zero
  Do i = 1, nev
     Call op(nrow, 1, evec(1:nrow, i), lde, avec, nrow)
     Call dgemv('T', nrow, i, one, evec, lde, avec, 1, zero,&
          & wrk2((i-1)*nev+1), 1)
  End Do
  Call trl_g_sum(info%mpicom, nsqr, wrk2(1:nsqr), uau)
  Do i = 2, nev
     wrk2(i:(i-1)*nev:nev) = wrk2((i-1)*nev+1:(i-1)*nev+i-1)
  End Do
  ! compute solution of R^T H_2 = H_1
  Call dtrtrs('U', 'T', 'N', nev, nev, rvv, nev, wrk2, nev, ierr)
  If (ierr .Ne. 0) Then
     info%stat = -235
     Goto 888
  End If
  ! compute solution of R^T H = H_2^T
  Do i = 1, nev
     uau(i:nsqr:nev) = wrk2((i-1)*nev+1:i*nev)
  End Do
  Call dtrtrs('U', 'T', 'N', nev, nev, rvv, nev, uau, nev, ierr)
  If (ierr .Ne. 0) Then
     info%stat = -236
     Goto 888
  End If
  ! solve the small symmetric eigenvalue problem
  Call dsyev('V', 'U', nev, uau, nev, eres, wrk2, nsqr, ierr)
  If (ierr .Ne. 0) Then
     info%stat = -237
     Goto 888
  End If
  ! solve R Y = Y to prepare for multiplying with V
  Call dtrtrs('U', 'N', 'N', nev, nev, rvv, nev, uau, nev, ierr)
  If (ierr .Ne. 0) Then
     info%stat = -238
     Goto 888
  End If
  ! call trl_ritz_vector to do the final multiplication
  If (lwrk .Ge. 3*nsqr) Then
     wrk2 => wrk(nsqr+1:lwrk)
  Else If (lwrk .Ge. nsqr+nsqr) Then
     wrk2 => wrk
  Else
     wrk2 => rvv
  End If
  i = Size(wrk2)
  Call trl_ritz_vectors(nrow, 0, nev, uau, nev,&
       & evec, lde, nev, avec, nrow, 0, wrk2, i)
  ! compute the residual norms
  Do i = 1, nev
     Call op(nrow, 1, evec(1:nrow,i), lde, avec, nrow)
     avec = avec - eres(i)*evec(1:nrow,i)
     eres(nev+i) = Dot_product(avec, avec)
  End Do
  Call trl_g_sum(info%mpicom, nev, eres(nev+1:nev+nev), avec)
  Do i = nev+1, nev+nev
     If (eres(i) .Gt. 0.0D0) Then
        eres(i) = Sqrt(eres(i))
     Else
        eres(i) = -Huge(eres(1))
     End If
  End Do
  If (info%lohi .Lt. 0) Then
     Do i = nev, nev+nev-2
        eres(i) = eres(i+1)
     End Do
  Else If (info%lohi .Gt. 0) Then
     Do i = 1, nev-1
        eres(i) = eres(i+1)
        evec(:,i) = evec(:,i+1)
     End Do
     Do i = nev, nev+nev-2
        eres(i) = eres(i+2)
     End Do
  End If
888 If (iuau .Eq. 0) Deallocate(uau)
  If (irvv .Eq. 0) Deallocate(rvv)
  If ((.Not. Present(base)) .And. mev.Le.nev) Deallocate(avec)
End Subroutine trl_ritz_projection
!!!
! a subroutine to compute Rayleigh quotients --
! Given a set of Ritz vectors and Ritz values, normalize the Ritz
! vectors and compute their Rayleigh quotients to replace the Ritz
! values.
!
! Arguments:
! op   -- the matrix-vector multiplication routine
! info -- the data structure that stores infomation for TRLAN
! evec -- the array to store the portion of eigenvectors on this PE
! base -- the workspace used to store results of MATVEC
! eres -- store new Ritz values and new Rresidual norms
!         if there are NEV Ritz pairs, eres(1:NEV) stores the new
!         Rayleigh quotient and eres(nev+1:nev+nev) stores the
!         new residual norms.
!!!
Subroutine trl_rayleigh_quotients(op, info, evec, eres, base)
  Use trl_info
  Implicit None
  Type(TRL_INFO_T) :: info
  Real(8), Dimension(:) :: eres
  Real(8), Dimension(:,:) :: evec
  Real(8), Dimension(:), Target, Optional :: base
  Interface
     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
  Integer :: i, nev, nrow
  Real(8), Dimension(4) :: wrk
  Real(8), Dimension(:), Pointer :: avec
  nev = Size(evec, 2)
  nrow = info%nloc
  If (nev .Le. 0) Return
  If (Present(base)) Then
     avec => base(1:nrow)
  Else
     Allocate(avec(nrow))
  End If
  avec = 0
  ! loop through each vector to normalize the vector, compute Rayleigh
  ! quotient and compute residual norm of the new Ritz pairs
  Do i = 1, nev
     wrk(1) = Dot_product(evec(1:nrow,i), evec(1:nrow,i))
     Call op(nrow, 1, evec(1:nrow,i), nrow, avec, nrow)
     wrk(2) = Dot_product(evec(1:nrow,i), avec)
     Call trl_g_sum(info%mpicom, 2, wrk(1:2), wrk(3:4))
     info%matvec = info%matvec + 1
     info%flop = info%flop + 4*nrow
     If (wrk(1) .Gt. 0.0D0) Then
        eres(i) = wrk(2)/wrk(1)
        avec = avec - eres(i)*evec(1:nrow,i)
        wrk(2) = Dot_product(avec, avec)
        Call trl_g_sum(info%mpicom, 1, wrk(2), wrk(3))
        wrk(1) = 1.0D0 / Sqrt(wrk(1))
        eres(nev+i) = wrk(1) * Sqrt(wrk(2))
        evec(1:nrow,i) = wrk(1) * evec(1:nrow,i)
        info%flop = info%flop + 6*nrow + 3
     Else
        eres(i) = -Huge(wrk(1))
        eres(nev+i) = -Huge(wrk(1))
     End If
  End Do
  If (.Not. Present(base)) Deallocate(avec)
End Subroutine trl_rayleigh_quotients