fstr_EIG_lanczos_util.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_lanczos_util
 
contains
 
  subroutine lanczos_set_initial_value(hecMESH, hecMAT, fstrEIG, eigvec, p, q, beta)
    use m_fstr
    use hecmw_util
    implicit none
    type(hecmwst_local_mesh) :: hecMESH
    type(hecmwst_matrix)     :: hecMAT
    type(fstr_eigen)         :: fstrEIG
    integer(kind=kint) :: N, NP, NDOF, NDOF2, NNDOF, NPNDOF
    integer(kind=kint) :: i, j
    real(kind=kreal)   :: eigvec(:, :), p(:), beta, chk, sigma
    real(kind=kreal), allocatable :: temp(:)
    real(kind=kreal), pointer     :: q(:), mass(:), filter(:)
 
    n      = hecmat%N
    np     = hecmat%NP
    ndof   = hecmesh%n_dof
    ndof2  = ndof*ndof
    nndof  = n *ndof
    npndof = np*ndof
 
    sigma  = fstreig%sigma
    mass   => fstreig%mass
    filter => fstreig%filter
 
    allocate(temp(nndof))
    temp = 0.0d0
 
    if(fstreig%is_free)then
      do i = 1, np
        do j = 1, ndof
          hecmat%D(ndof2*(i-1) + (ndof+1)*(j-1) + 1) = hecmat%D(ndof2*(i-1) + (ndof+1)*(j-1) + 1) + sigma * mass(ndof*(i-1) + j)
        enddo
      enddo
    endif
 
    call urand1(nndof, temp, hecmesh%my_rank)
 
    do i = 1, nndof
      temp(i) = temp(i) * filter(i)
    enddo
 
    do i = 1, nndof
      eigvec(i,1) = mass(i) * temp(i)
    enddo
 
    chk = 0.0d0
    do i = 1, nndof
      chk = chk + temp(i) * eigvec(i,1)
    enddo
    call hecmw_allreduce_r1(hecmesh, chk, hecmw_sum)
    beta = dsqrt(chk)
 
    if(beta == 0.0d0)then
      call hecmw_finalize()
      stop "Self-orthogonal"
    endif
 
    chk = 1.0d0/beta
    do i = 1, nndof
      q(i) = temp(i) * chk
    enddo
 
    do i = 1, nndof
      p(i) = mass(i) * q(i)
    enddo
  end subroutine lanczos_set_initial_value
 
  subroutine evsort(EIG, NEW, NEIG)
    use hecmw
    implicit none
    integer(kind=kint) :: i, j, n, ip, minloc, NEIG, IBAF, NEW(NEIG)
    real(kind=kreal) :: emin, eig(neig)
 
    do i = 1, neig
      new(i) = i
    enddo
 
    n = neig-1
    do i = 1, n
      minloc = i
      emin = dabs(eig(new(i)))
      ip = i+1
      do j=ip,neig
        if(dabs(eig(new(j))).LT.emin)then
          minloc = j
          emin = dabs(eig(new(j)))
        endif
      enddo
      ibaf = new(i)
      new(i) = new(minloc)
      new(minloc) = ibaf
    enddo
  end subroutine evsort
 
  subroutine urand1(N, X, SHIFT)
    use hecmw
    implicit none
    real(kind=kreal) :: x(n), invm
    integer(kind=kint), parameter :: MM = 1664501
    integer(kind=kint), parameter :: LAMBDA = 1229
    integer(kind=kint), parameter :: MU = 351750
    integer(kind=kint) :: i, N, IR, SHIFT
 
    ir = 0
    invm = 1.0d0 / mm
    do i = 1, shift
      ir = mod( lambda * ir + mu, mm)
    enddo
    do i = shift+1, shift+n
      ir = mod( lambda * ir + mu, mm)
      x(i-shift) = invm * ir
    enddo
  end subroutine urand1
 
end module m_fstr_eig_lanczos_util