fstr_EIG_lanczos.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
contains
 
  subroutine fstr_solve_lanczos(hecMESH, hecMAT, fstrSOLID, fstrEIG)
    use m_fstr
    use hecmw_util
    use m_eigen_lib
    use m_fstr_eig_lanczos_util
    use m_fstr_eig_tridiag
 
    implicit none
 
    type(hecmwst_local_mesh) :: hecmesh
    type(hecmwst_matrix)     :: hecMAT
    type(fstr_solid)         :: fstrSOLID
    type(fstr_eigen)         :: fstrEIG
    type(fstr_tri_diag)      :: Tri
    type(fstr_eigen_vec), pointer :: Q(:)
    integer(kind=kint) :: N, NP, NDOF, NNDOF, NPNDOF
    integer(kind=kint) :: iter, maxiter, nget, ierr
    integer(kind=kint) :: i, j, k, in, jn, kn, ik, it
    integer(kind=kint) :: ig, ig0, is0, ie0, its0, ite0
    real(kind=kreal)   :: t1, t2, tolerance
    real(kind=kreal)   :: alpha, beta, beta0
    real(kind=kreal), allocatable :: s(:), t(:), p(:)
    logical :: is_converge
 
    n      = hecmat%N
    np     = hecmat%NP
    ndof   = hecmesh%n_dof
    nndof  = n *ndof
    npndof = np*ndof
 
    allocate(fstreig%filter(npndof))
    fstreig%filter = 1.0d0
    fstreig%sigma = 0.01d0
 
    jn = 0
    do ig0 = 1, fstrsolid%BOUNDARY_ngrp_tot
      ig   = fstrsolid%BOUNDARY_ngrp_ID(ig0)
      is0  = hecmesh%node_group%grp_index(ig-1) + 1
      ie0  = hecmesh%node_group%grp_index(ig  )
      it   = fstrsolid%BOUNDARY_ngrp_type(ig0)
      its0 = (it - mod(it,10))/10
      ite0 = mod(it,10)
 
      do ik = is0, ie0
        in = hecmesh%node_group%grp_item(ik)
        if(ndof < ite0) ite0 = ndof
        do i = its0, ite0
          jn = jn + 1
          fstreig%filter((in-1)*ndof+i) = 0.0d0
        enddo
      enddo
    enddo
 
    do ig0 = 1, fstrsolid%SPRING_ngrp_tot
      ig = fstrsolid%SPRING_ngrp_ID(ig0)
      is0 = hecmesh%node_group%grp_index(ig-1) + 1
      ie0 = hecmesh%node_group%grp_index(ig  )
      do ik = is0, ie0
        jn = jn + 1
      enddo
    enddo
 
    call hecmw_allreduce_i1(hecmesh, jn, hecmw_sum)
    if(jn == 0)then
      fstreig%is_free = .true.
      if(myrank == 0)then
        write(*,"(a,1pe12.4)") '** free modal analysis: shift factor =', fstreig%sigma
      endif
    endif
 
    call hecmw_update_r(hecmesh, fstreig%filter, np, ndof)
 
    in = 0
    do i = 1, nndof
      if(fstreig%filter(i) == 1.0d0) in = in + 1
    enddo
    call hecmw_allreduce_i1(hecmesh, in, hecmw_sum)
 
    fstreig%maxiter = fstreig%maxiter + 1
    if(in < fstreig%maxiter)then
      if(myrank == 0)then
        write(imsg,*) '** changed maxiter to system matrix size.'
      endif
      fstreig%maxiter = in
    endif
 
    if(in < fstreig%nget)then
      fstreig%nget = in
    endif
 
    maxiter   = fstreig%maxiter
 
    allocate(q(0:maxiter))
    allocate(q(0)%q(npndof))
    allocate(q(1)%q(npndof))
    allocate(fstreig%eigval(maxiter))
    allocate(fstreig%eigvec(npndof, maxiter))
    allocate(tri%alpha(maxiter))
    allocate(tri%beta(maxiter))
    allocate(t(npndof))
    allocate(s(npndof))
    allocate(p(npndof))
 
    fstreig%eigval = 0.0d0
    fstreig%eigvec = 0.0d0
    t      = 0.0d0
    p      = 0.0d0
    s      = 0.0d0
    q(0)%q = 0.0d0
    q(1)%q = 0.0d0
    tri%alpha = 0.0d0
    tri%beta  = 0.0d0
    hecmat%X  = 0.0d0
 
    call lanczos_set_initial_value(hecmesh, hecmat, fstreig, fstreig%eigvec, p, q(1)%q, tri%beta(1))
 
    hecmat%Iarray(98) = 1 !Assembly complete
    hecmat%Iarray(97) = 1 !Need numerical factorization
 
    if(myrank == 0)then
      write(imsg,*)
      write(imsg,*) ' *****   STAGE Begin Lanczos loop     **'
    endif
 
    do iter = 1, maxiter-1
      do i = 1, npndof
        hecmat%B(i) = p(i)
      enddo
 
      call solve_lineq(hecmesh, hecmat)
 
      allocate(q(iter+1)%q(npndof))
 
      do i = 1, npndof
        t(i) = hecmat%X(i) * fstreig%filter(i)
      enddo
 
      do i = 1, npndof
        t(i) = t(i) - tri%beta(iter) * q(iter-1)%q(i)
      enddo
 
      alpha = 0.0d0
      do i = 1, nndof
        alpha = alpha + p(i) * t(i)
      enddo
      call hecmw_allreduce_r1(hecmesh, alpha, hecmw_sum)
      tri%alpha(iter) = alpha
 
      do i = 1, npndof
        t(i) = t(i) - tri%alpha(iter) * q(iter)%q(i)
      enddo
 
      s = 0.0d0
 
      do i = 1, npndof
        s(i) = fstreig%mass(i) * t(i)
      enddo
 
      do j = 0, iter
        t1 = 0.0d0
        do i = 1, nndof
          t1 = t1 + q(j)%q(i) * s(i)
        enddo
        call hecmw_allreduce_r1(hecmesh, t1, hecmw_sum)
        do i = 1, npndof
          t(i) = t(i) - t1 * q(j)%q(i)
        enddo
      enddo
 
      do i = 1, npndof
        s(i) = fstreig%mass(i) * t(i)
      enddo
 
      beta = 0.0d0
      do i = 1, nndof
        beta = beta + s(i) * t(i)
      enddo
      call hecmw_allreduce_r1(hecmesh, beta, hecmw_sum)
      tri%beta(iter+1) = dsqrt(beta)
 
      beta = 1.0d0/tri%beta(iter+1)
      do i = 1, npndof
        p(i)           = s(i) * beta
        q(iter+1)%q(i) = t(i) * beta
      enddo
 
      fstreig%iter = iter
      if(iter == 1) beta0 = tri%beta(iter+1)
 
      call tridiag(hecmesh, hecmat, fstreig, q, tri, iter, is_converge)
 
      if(is_converge) exit
    enddo
 
    do i = 0, iter
      if(associated(q(i)%q)) deallocate(q(i)%q)
    enddo
    deallocate(tri%alpha)
    deallocate(tri%beta)
    deallocate(t)
    deallocate(s)
    deallocate(p)
 
    t2 = hecmw_wtime()
 
    if(myrank == 0)then
      write(imsg,*)
      write(imsg,*) ' *     STAGE Output and postprocessing    **'
      write(idbg,'(a,f10.2)') 'Lanczos loop (sec) :', t2 - t1
    endif
 
  end subroutine fstr_solve_lanczos
 
end module m_fstr_eig_lanczos