trlcore.f90

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

           res(i) = res(j)
           lambda(i) = lambda(j)
           j = j + 1
        End Do
     End If
  End If
End Subroutine trl_sort_eig
!!!
! generating eigenvectors of the projected eigenvalue problem
! corresponding to the given Ritz values
! using LAPACK routine DSTEIN (inverse iterations)
!
! workspace size:
! iwrk: dimension(4*nd)
! wrk:  lwrk >= 5*nd
!!!
Subroutine trl_get_tvec(nd, alpha, beta, irot, nrot, rot, nlam,&
     & lambda, yy, iwrk, wrk, lwrk, ierr)
  Implicit None
  Integer, Intent(in) :: nd, nrot, nlam, irot, lwrk
  Integer :: ierr, iwrk(nd,4)
  Real(8), Intent(in) :: alpha(nd), beta(nd), lambda(nlam), rot(nrot*nrot)
  Real(8) :: wrk(lwrk), yy(nd*nlam)
  !
  ! local variables
  Character, Parameter :: notrans='N'
  Real(8), Parameter :: zero=0.0D0, one=1.0D0
  Integer :: i, j, k, ncol, ii, ioff
  !
  ! conventional external subprograms
  External dstein, dgemm, dgemv
  If (nlam .Le. 0) Return
  If (lwrk.Ge.5*nd) Then
     ierr = 0
  Else
     ierr = -131
     Return
  End If
  !
  ! set up IBLOCK and ISPLIT for calling dstein
  iwrk(1:nd,1) = 1
  iwrk(1:nd,2) = nd
  Call dstein(nd, alpha, beta, nlam, lambda, iwrk(1,1), iwrk(1,2), yy,&
       & nd, wrk, iwrk(1,3), iwrk(1,4), ierr)
  If (ierr .Ne. 0) Then
     Print *, 'TRL_GET_TVEC: dstein failed with error code ', ierr
     ierr = -132
     Return
  End If
  !
  ! apply the rotations to the IROT+1:IROT+NROT rows of YY
  ! generates results 'NCOL' columns at a time
  If (nrot.Gt.1) Then
     ncol = lwrk / nrot
     Do i = 1, nlam, ncol
        j = Min(nlam, i+ncol-1)
        k = j - i + 1
        If (k .Gt. 1) Then
           Call dgemm(notrans, notrans, nrot, k, nrot, one, rot, nrot,&
                & yy((i-1)*nd+irot+1), nd, zero, wrk, nrot)
           Do ii = i-1, j-1
              ioff = (ii-i+1)*nrot
              yy(ii*nd+irot+1:ii*nd+irot+nrot) = wrk(ioff+1:ioff+nrot)
           End Do
        Else
           Call dgemv(notrans, nrot, nrot, one, rot, nrot,&
                & yy((i-1)*nd+irot+1), 1, zero, wrk, 1)
           yy((i-1)*nd+irot+1:(i-1)*nd+irot+nrot) = wrk(1:nrot)
        End If
     End Do
  End If
End Subroutine trl_get_tvec
!!!
! compute all eigenvalues and eigenvectors of the projected matrix
! use LAPACK routine DSYEV
! The eigenvectors corresponding to lambda(1:nlam) are placed at the
! first nlam*nd locations of yy on exit.
!
! On entry, Alpha and Beta defines the arrayhead plus tridiagonal
! matrix.
! On return, Alpha will store the Ritz values in the same order as
! lambda.
!
! lwrk >= 3*nd
!!!
Subroutine trl_get_tvec_a(nd, kept, alpha, beta, nlam, lambda,&
     & yy, wrk, lwrk, iwrk, ierr)
  Implicit None
  Integer, Intent(in) :: kept, nd, nlam, lwrk
  Integer :: ierr, iwrk(nd)
  Real(8), Intent(in) :: lambda(nd)
  Real(8) :: alpha(nd), beta(nd), yy(nd*nd), wrk(lwrk)
  !
  ! local variables
  Integer :: i, j, i2, j2, ii
  Real(8) :: tmp
  External dsyev
  !
  ! fill yy with the projection matrix, then call DSYEV
  If (nlam .Le. 0) Return
  If (lwrk.Ge.nd+nd+nd) Then
     ierr = 0
  Else
     ierr = -141
     Return
  End If
  yy = 0.0D0
  yy(1:nd*nd:nd+1) = alpha(1:nd)
  If (kept.Gt.0) yy(kept*nd+1:kept*nd+kept) = beta(1:kept)
  yy((kept+1)*(nd+1):nd*nd:nd+1) = beta(kept+1:nd-1)
  Call dsyev('V', 'U', nd, yy, nd, alpha, wrk, lwrk, ierr)
  If (ierr .Ne. 0) Then
     ierr = -142
     Return
  End If
  If (nlam .Ge. nd) Return
  !
  ! reorder the eigenvectors
  ! both lambda(1:kept) and alpha are in ascending order
  !
  tmp = Max(alpha(nd)-alpha(1), Abs(alpha(nd)), Abs(alpha(1)))
  tmp = Epsilon(tmp)*tmp*nd
  j = 1
  i = 1
  Do While (i .Le. nlam)
     ! move j so that alpha(j) is within tmp distance away
     ii = j
     j = nd
     Do While (ii .Le. nd)
        If (alpha(ii) .Lt. lambda(i)-tmp) Then
           ii = ii + 1
        Else
           j = ii
           ii = nd+1
        End If
     End Do
     If (alpha(j) .Gt. lambda(i)+tmp) Then
        ierr = -143
        Return
     End If
     ! identify the group size in lambda
     ii = i + 1
     i2 = nlam
     Do While (ii .Le. nlam)
        If (lambda(ii) .Le. lambda(i)+tmp) Then
           ii = ii + 1
        Else
           i2 = ii - 1
           ii = nd+1
        End If
     End Do
     ! identify the group size in alpha
     ii = j + 1
     j2 = nd
     Do While (ii .Le. nd)
        If (alpha(ii) .Le. lambda(i)+tmp) Then
           ii = ii + 1
        Else
           j2 = ii - 1
           ii = nd+1
        End If
     End Do
     ! assign the index values
     If (j2.Eq.j .And. i2.Eq.i) Then
        iwrk(i) = j
     Else If (j2-j .Eq. i2-i) Then
        iwrk(i:i2) = (/(ii, ii=j, j2)/)
     Else If (j2-j .Gt. i2-i) Then
        j2 = j + i2 - i
        iwrk(i:i2) = (/(ii, ii=j, j2)/)
     Else If (j2 .Lt. nd) Then
        i2 = i + j2 - j
        iwrk(i:i2) = (/(ii, ii=j, j2)/)
     Else
        ierr = -144
        Return
     End If
     i = i2 + 1
     j = j2 + 1
  End Do
  ! perform the actual copying
  Do i = 1, nlam
     j = iwrk(i)
     If (j.Gt.i) Then
        alpha(i) = alpha(j)
        yy((i-1)*nd+1:i*nd) = yy((j-1)*nd+1:j*nd)
     End If
  End Do
End Subroutine trl_get_tvec_a
!!!
! move the Ritz pairs with extremely small residual norms to the front
! of the arrays so that locking can be performed cleanly
!!!
Subroutine trl_set_locking(jnd, nlam, lambda, res, yy, anrm, locked)
  Implicit None
  Integer, Intent(in) :: jnd, nlam
  Integer :: locked
  Real(8), Intent(in) :: anrm
  Real(8) :: lambda(nlam), res(nlam), yy(jnd*nlam)
  !
  Real(8), Parameter :: zero = 0.0D0
  Integer :: i, j, ii, ioff
  Real(8) :: tmp, eps, small
  Logical :: ti, tj
  small(tmp, eps) = eps*Max(Abs(tmp), eps*anrm)
  eps = Epsilon(zero)
  i = 1
  j = nlam
  ti = (Abs(res(i)) .Le. small(lambda(i),eps))
  tj = (Abs(res(j)) .Le. small(lambda(j),eps))
  Do While (i .Lt. j)
     If (ti) Then
        res(i) = zero
        i = i + 1
        If (i.Le.j) Then
           ti = (Abs(res(i)) .Le. small(lambda(i),eps))
        Else
           ti = .False.
        End If
     Else
        If (tj) Then
           tmp = lambda(i)
           lambda(i) = lambda(j)
           lambda(j) = tmp
           res(j) = res(i)
           res(i) = zero
           ioff = (j-i)*jnd
           Do ii = i*jnd-jnd+1, i*jnd
              tmp = yy(ii)
              yy(ii) = yy(ii+ioff)
              yy(ii+ioff) = tmp
           End Do
           i = i + 1
           If (i.Le.j) Then
              ti = (Abs(res(i)) .Le. small(lambda(i),eps))
           Else
              ti = .False.
           End If
        End If
        j = j - 1
        If (j.Gt.i) Then
           tj = (Abs(res(j)) .Le. small(lambda(j),eps))
        Else
           tj = .False.
        End If
     End If
  End Do
  If (ti) Then
     locked = i
  Else
     locked = i - 1
  End If
End Subroutine trl_set_locking
!!!
! compute the Ritz vectors from the basis vectors and the eigenvectors
! of the projected system
! the basis vectors may be stored in two separete arrays
! the result need to be stored back in them
!
! lwrk should be no less than ny (lwrk>=ny) ***NOT checked inside***
!!!
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)
  !
  ! local variables
  !
  Character, Parameter :: notrans='N'
  Real(8), Parameter :: zero=0.0D0, one=1.0D0
  Integer :: i,j,k,stride,ii,jl1,jl2,il1,il2,kv1
  !
  ! conventional external procedures
  External dgemm, dgemv
  !
  If (lck .Le. m1) Then
     il1 = lck + 1
     jl1 = m1 - lck
     il2 = 1
     jl2 = m2
  Else
     il1 = m1+1
     jl1 = 0
     il2 = lck - m1 + 1
     jl2 = m1 + m2 - lck
  End If
  If (jl1.Eq.0 .And. jl2.Eq.0) Return
  kv1 = Min(m1-il1+1, ny)
  wrk = zero
  If (ny .Gt. 1) Then
     stride = lwrk / ny
     Do i = 1, nrow, stride
        j = Min(nrow, i+stride-1)
        k = j - i + 1
        If (jl1 .Gt. 0) Then
           Call dgemm(notrans, notrans, k, ny, jl1, one, vec1(i,il1),&
                & ld1, yy, ldy, zero, wrk, k)
        Else
           wrk = zero
        End If
        If (jl2 .Gt. 0) Call dgemm(notrans, notrans, k, ny, jl2, one,&
             & vec2(i,il2), ld2, yy(jl1+1,1), ldy, one, wrk, k)
        Do ii = 0, kv1-1
           vec1(i:j,ii+il1) = wrk(ii*k+1:ii*k+k)
        End Do
        Do ii = 0, ny-kv1-1
           vec2(i:j,ii+il2) = wrk((kv1+ii)*k+1:(kv1+ii)*k+k)
        End Do
     End Do
  Else If (ny.Eq.1) Then
     stride = lwrk
     Do i = 1, nrow, stride
        j = Min(nrow, i+stride-1)
        k = j - i + 1
        If (jl1 .Gt. 0) Then
           Call dgemv(notrans, k, jl1, one, vec1(i,il1), ld1, yy, 1,&
                & zero, wrk, 1)
           If (jl2 .Gt. 0) Call dgemv(notrans, k, jl2, one,&
                & vec2(i,il2), ld2, yy(jl1+1,1), 1, one, wrk, 1)
        Else
           Call dgemv(notrans, k, jl2, one, vec2(i,il2),&
                & ld2, yy(jl1+1,1), 1, zero, wrk, 1)
        End If
        If (kv1 .Gt. 0) Then
           vec1(i:j,il1) = wrk(1:k)
        Else
           vec2(i:j,il2) = wrk(1:k)
        End If
     End Do
  End If
End Subroutine trl_ritz_vectors
!!!
! full Gram-Schmidt routine -- orthogonalize a new vector against
! all existing vectors
!
! On entry:
! rnrm and alpha are expected to have the correct values.
!!!
Subroutine trl_cgs(mpicom, myid, nrow, v1, ld1, m1, v2, ld2, m2, rr,&
     & rnrm, alpha, north, nrand, wrk, ierr)
  Implicit None
  Integer, Intent(in) :: mpicom, myid, nrow, ld1, m1, ld2, m2
  Integer, Intent(inout) :: north, nrand, ierr
  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)
  Interface
     ! global dot-product routine, calls dgemv
     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
     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
  !
  Character, Parameter :: notrans='N'
  Integer, Parameter :: maxorth = 5
  Real(8), Parameter :: one=1.0D0, zero=0.0D0
  Integer :: i, j, nold, irnd, cnt
  Real(8) :: tmp
  Logical :: sane
  External dgemv
  !
  nold = m1+m2
  If (ld1.Ge.nrow .And. (ld2.Ge.nrow.Or.m2.Le.0)) Then
     ierr = 0
  Else
     ierr = -201
     Return
  End If
  irnd = 0
  If (nold .Gt. 0) Then
     ! repeated performing Gram-Schmidt procedure to ensure
     ! ortho