fstr_EIG_tridiag.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
module m_fstr_eig_tridiag
  use hecmw
 
  implicit none
 
  public
 
  type fstr_eigen_vec
    real(kind=kreal), pointer :: q(:) => null()
  end type fstr_eigen_vec
 
  type fstr_tri_diag
    real(kind=kreal), allocatable :: alpha(:)
    real(kind=kreal), allocatable :: beta(:)
  end type fstr_tri_diag
 
contains
 
  subroutine tridiag(hecMESH, hecMAT, fstrEIG, Q, Tri, iter, is_converge)
    use hecmw
    use m_fstr
    use m_fstr_eig_lanczos_util
    implicit none
    type(hecmwst_local_mesh) :: hecmesh
    type(hecmwst_matrix)     :: hecMAT
    type(fstr_eigen)         :: fstrEIG
    type(fstr_tri_diag)      :: Tri
    type(fstr_eigen_vec), pointer :: Q(:)
 
    integer(kind=kint), intent(in) :: iter
    integer(kind=kint) :: N, NP, NDOF, NNDOF, NPNDOF
    integer(kind=kint) :: i, j, k, in, jn, kn, nget
    integer(kind=kint) :: iter2, ierr, maxiter
    real(kind=kreal)   :: resid, chk, vmax, sigma, tolerance
    real(kind=kreal), allocatable :: alpha(:), beta(:), temp(:)
    real(kind=kreal), allocatable :: l(:,:)
 
    integer(kind=kint), allocatable :: iparm(:)
    real(kind=kreal), pointer       :: eigvec(:,:)
    real(kind=kreal), pointer       :: eigval(:)
    logical :: is_converge
 
    n      = hecmat%N
    np     = hecmat%NP
    ndof   = hecmesh%n_dof
    nndof  = n *ndof
    npndof = np*ndof
    eigval => fstreig%eigval
    eigvec => fstreig%eigvec
    nget      = fstreig%nget
    maxiter   = fstreig%maxiter
    tolerance = fstreig%tolerance
 
    allocate( iparm(maxiter) )
    allocate( temp(maxiter)  )
    allocate( alpha(iter)    )
    allocate( beta(iter)     )
    allocate( l(iter, iter)  )
 
    do j=1, iter
      do i = 1,iter
        l(i,j)  = 0.0d0
      enddo
    enddo
 
    do i=1, iter
      alpha(i) = tri%alpha(i)
      l(i,i)   = 1.0d0
    enddo
 
    beta(1) = 0.0d0
    do i=2, iter
      beta(i) = tri%beta(i)
    enddo
 
    call ql_decomposition(iter, iter, alpha, beta, l, ierr)
 
    sigma = 0.0d0
    if(fstreig%is_free) sigma = fstreig%sigma
    do i = 1, iter
      if(alpha(i) /= 0.0d0)then
        eigval(i) = 1.0d0/alpha(i) - sigma
      else
        ! A zero Ritz value stands for an infinite eigenvalue. evsort orders by
        ! magnitude, so the sentinel keeps such a mode out of the requested set.
        eigval(i) = huge(0.0d0)
      endif
    enddo
 
    ! The lowest modes are the largest Ritz values, so the convergence check has to
    ! pick them through iparm rather than take the leading entries of alpha.
    call evsort(eigval, iparm, iter)
 
    is_converge = .true.
    chk = 0.0d0
    ! Two extra modes absorb the reordering of modes sitting on the boundary of the requested set.
    do i = 1, min(nget+2, iter)
      in = iparm(i)
      if (dabs(alpha(in)) > 0.0d0) then
        ! Divide by the eigenvalue of the mode itself.
        resid = dabs(tri%beta(iter+1)*l(iter,in))/dabs(alpha(in))
        chk = max(chk, resid)
        if(tolerance < resid) is_converge = .false.
      else
        is_converge = .false.
      endif
    enddo
    if(myrank == 0) write(*,"(i8,1pe12.5)")iter, chk
 
    if(iter < nget) is_converge = .false.
 
    if(iter == maxiter-1 .and. .not. is_converge)then
      if(myrank == 0)then
        write(*,*)    '### WARNING: eigen analysis stopped at maxiter without convergence.'
        write(ilog,*) '### WARNING: eigen analysis stopped at maxiter without convergence.'
      endif
      ! Only iter Ritz pairs exist. Reporting more would print entries of eigval and
      ! eigvec that no iteration has ever written.
      if(iter < nget) fstreig%nget = iter
      is_converge = .true.
    endif
 
    if(is_converge)then
      temp = eigval
 
      eigvec = 0.0d0
      do k=1, iter
        in = iparm(k)
        eigval(k) = temp(in)
        do j=1, iter
          do i=1, npndof
            eigvec(i, k) = eigvec(i, k) + q(j)%q(i) * l(j, in)
          enddo
        enddo
      enddo
 
      do j=1, iter
        chk  = 0.0d0
        vmax = 0.0d0
        do i = 1, nndof
          chk  = max(chk,  dabs(eigvec(i,j)))
          vmax = max(vmax, eigvec(i,j))
        enddo
        call hecmw_allreduce_r1(hecmesh, chk, hecmw_max)
        call hecmw_allreduce_r1(hecmesh, vmax, hecmw_max)
        if(chk /= 0.0d0)then
          ! The largest component becomes +1. A tie between +a and -a goes to the
          ! positive one so that serial and parallel runs pick the same sign.
          if(vmax /= chk) chk = -chk
          chk = 1.0d0/chk
          do i = 1, nndof
            eigvec(i,j) = eigvec(i,j) * chk
          enddo
        endif
      enddo
    endif
 
    deallocate(iparm)
    deallocate(temp)
    deallocate(alpha)
    deallocate(beta)
    deallocate(l)
  end subroutine tridiag
 
  !======================================================================!
  !                       Description                                    !
  !======================================================================!
  !
  !on input
  !
  !nm must be set to the row dimension of two-dimensional
  !array parameters as declared in the calling program
  !dimension statement.
  !
  !is the order of the matrix.
  !
  !contains the diagonal elements of the input matrix.
  !
  !contains the subdiagonal elements of the input matrix
  !in its last n-1 positions.  e(1) is arbitrary.
  !
  !contains the transformation matrix produced in the
  !reduction by  tred2, if performed.  if the eigenvectors
  !of the tridiagonal matrix are desired, z must contain
  !the identity matrix.
  !
  !on output
  !
  !d
  !contains the eigenvalues in the order produced by the ql sweeps,
  !which is not sorted.  ordering is left to the caller.  if an
  !error exit is made, the eigenvalues are correct for indices
  !1,2,...,ierror-1.
  !-----------------------------------------------------
  !e
  !has been destroyed.
  !
  !z
  !contains orthonormal eigenvectors of the symmetric
  !tridiagonal (or full) matrix.  if an error exit is made,
  !z contains the eigenvectors associated with the stored
  !eigenvalues.
  !-----------------------------------------------------
  !ierror is set to
  !  zero       for normal return,
  !  j          if the j-th eigenvalue has not been
  !             determined after 30 iterations.
  !
  !calls a2b2 for  dsqrt(a*a + b*b) .
  !=======================================================================
 
  !call QL_decomposition(iter, iter, alpha, beta, L, ierr)
  subroutine ql_decomposition(nm, n, d, e, z, ierror)
    use hecmw
    implicit none
    integer(kind=kint) :: i, j, k, l, m, n, ii, l1, l2, nm, mml, ierror
    real(kind=kreal)   :: d(n), e(n), z(nm, n)
    real(kind=kreal)   :: c, c2, c3, dl1, el1, f, g, h, p, r, s, s2, tst1, tst2
 
    ierror = 0
    if (n .eq. 1) go to 1001
 
    do i = 2, n
      e(i-1) = e(i)
    enddo
 
    f = 0.0d0
    tst1 = 0.0d0
    e(n) = 0.0d0
 
    do 240 l = 1, n
      j = 0
      h = dabs(d(l)) + dabs(e(l))
      if (tst1 .lt. h) tst1 = h
      !     .......... look for small sub-diagonal element ..........
      bb:do m = l, n
        tst2 = tst1 + dabs(e(m))
        if (tst2 .eq. tst1) exit bb
        !     .......... e(n) is always zero, so there is no exit
        !                through the bottom of the loop ..........
      enddo bb
 
      if (m .eq. l)  go to 220
 
      130    if (j .eq. 30) go to 1000
      j = j + 1
      !     .......... form shift ..........
      l1 = l + 1
      l2 = l1 + 1
      g = d(l)
      p = (d(l1) - g) / (2.0d0 * e(l))
      r = a2b2(p,1.0d0)
      d(l) = e(l) / (p + dsign(r,p))
      d(l1) = e(l) * (p + dsign(r,p))
      dl1 = d(l1)
      h = g - d(l)
      if (l2 .gt. n) go to 145
 
      do i = l2, n
        d(i) = d(i) - h
      enddo
 
      145    f = f + h
      !     .......... ql transformation ..........
      p = d(m)
      c = 1.0d0
      c2 = c
      el1 = e(l1)
      s = 0.0d0
      s2 = 0.0d0
      c3 = 0.0d0
      mml = m - l
      !     .......... for i=m-1 step -1 until l do -- ..........
      do ii = 1, mml
        c3 = c2
        c2 = c
        s2 = s
        i = m - ii
        g = c * e(i)
        h = c * p
        r = a2b2(p,e(i))
        e(i+1) = s * r
        s = e(i) / r
        c = p / r
        p = c * d(i) - s * g
        d(i+1) = h + s * (c * g + s * d(i))
        !     .......... form vector ..........
        do k = 1, n
          h = z(k,i+1)
          z(k,i+1) = s * z(k,i) + c * h
          z(k,i) = c * z(k,i) - s * h
        enddo
      enddo
 
      p = -s * s2 * c3 * el1 * e(l) / dl1
      e(l) = s * p
      d(l) = c * p
      tst2 = tst1 + dabs(e(l))
      if (tst2 .gt. tst1) go to 130
      220    d(l) = d(l) + f
      240 continue
 
      go to 1001
      !     .......... set error -- no convergence to an
      !                eigenvalue after 30 iterations ..........
      1000 ierror = l
      1001 return
  end subroutine ql_decomposition
 
  function a2b2(a,b)
    use hecmw
    implicit none
 
    real(kind=kreal) :: a2b2
    real(kind=kreal) :: a, b
    real(kind=kreal) :: p, q, r, s, t, u
 
    p = dmax1(dabs(a), dabs(b))
    if (p /= 0.0d0) then
      r = (dmin1(dabs(a),dabs(b))/p) ** 2
      do
        t = 4.0d0 + r
        if (t == 4.0d0) exit
        s = r/t
        u = 1.0d0 + 2.0d0*s
        p = u * p
        q = s/u
        r = q * q * r
      end do
    end if
    a2b2 = p
    return
 
  end function a2b2
 
end module m_fstr_eig_tridiag