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, sigma, tolerance, max_eigval
    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)
 
    is_converge = .true.
    chk = 0.0d0
    max_eigval = maxval(alpha)
    do i = 1, min(nget, iter)
      resid = dabs(tri%beta(iter+1)*l(iter,i))/max_eigval
      chk = max(chk, resid)
      if(tolerance < resid) is_converge = .false.
    enddo
    if(myrank == 0) write(*,"(i8,1pe12.5)")iter, chk
 
    if(iter < nget) is_converge = .false.
    if(iter == maxiter-1) is_converge = .true.
 
    if(is_converge)then
      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
        endif
      enddo
 
      call evsort(eigval, iparm, iter)
 
      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 = maxval(eigvec(:,j))
        call hecmw_allreduce_r1(hecmesh, chk, hecmw_max)
        if(chk /= 0.0d0)then
          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 ascending order.  if an
  !error exit is made, the eigenvalues are correct but
  !unordered for indices 1,2,...,ierror-1.
  !-----------------------------------------------------
  !GP
  !du
  !contains the unordered eigenvalues.  if an
  !error exit is made, the eigenvalues are correct and
  !unordered 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.
  !-----------------------------------------------------
  !zu
  !contains unordered eigenvectors of the symm. tridiag. matrix
  !-----------------------------------------------------
  !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
 
      !     .......... order eigenvalues and eigenvectors ..........
      !GP: Get unordered eigenvalues and eigenvectors----------------
 
      do 300 ii = 2, n
        i = ii - 1
        k = i
        p = d(i)
 
        aa:do j = ii, n
          if (d(j) .ge. p) exit aa
          k = j
          p = d(j)
        enddo aa
 
        if (k .eq. i) go to 300
        d(k) = d(i)
        d(i) = p
 
        do j = 1, n
          p = z(j,i)
          z(j,i) = z(j,k)
          z(j,k) = p
        enddo
 
        300 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