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, hecEBC)
    use m_fstr
    use hecmw_util
    use m_eigen_lib
    use m_fstr_eig_lanczos_util
    use m_fstr_eig_tridiag
    use hecmw_ebc_defer
    use hecmw_solver_las
 
    implicit none
 
    type(hecmwst_local_mesh) :: hecMESH
    type(hecmwst_matrix)     :: hecMAT
    type(fstr_solid)         :: fstrSOLID
    type(fstr_eigen)         :: fstrEIG
    type(hecmwst_ebc)        :: hecEBC
    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
    integer(kind=kint) :: ig, ig0, is0, ie0
    real(kind=kreal)   :: t1, t2, tolerance
    real(kind=kreal)   :: alpha, beta, beta0, resid
    real(kind=kreal), allocatable :: s(:), t(:), p(:)
    integer(kind=kint), allocatable :: mark(:)
    character(len=HECMW_MSG_LEN) :: msg(2), noconv
    logical :: is_converge, unusable
    ! a linear solution whose relative residual reaches this carries an error of about one percent, and
    ! the eigenvalues computed from it inherit that error; falling short of the tolerance by less than
    ! this is common and leaves the eigenvalues usable
    real(kind=kreal), parameter :: resid_unusable = 1.0d-2
 
    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 i = 1, npndof
      if(hecebc%mark(i) == 0) cycle
      jn = jn + 1
      fstreig%filter(i) = 0.0d0
    enddo
 
    if(hecmw_mat_get_mpc_method(hecmat) == 3)then
      allocate(mark(npndof))
      call hecmw_mpc_mark_slave(hecmesh, hecmat, mark)
      do i = 1, npndof
        if(mark(i) == 1) fstreig%filter(i) = 0.0d0
      enddo
      deallocate(mark)
    endif
 
    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 + 1
    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)
 
      ! the Lanczos vectors are built from this solution, so a solver that fell short of its tolerance is
      ! tolerated only while the solution stays usable
      if(hecmw_mat_get_flag_diverged(hecmat) /= kno .or. hecmw_mat_get_flag_converged(hecmat) == kno)then
        resid = hecmw_rel_resid_l2(hecmesh, hecmat)
        unusable = (resid >= resid_unusable .or. resid /= resid)  ! the second test catches NaN
        write(noconv,'(a,i0,a,1pe12.5)') ' the linear solver did not converge at Lanczos iteration ', &
          & iter, '; relative residual =', resid
        if(unusable)then
          msg(1) = '### ERROR:'//trim(noconv)
          msg(2) = '   the eigenvalues cannot be computed from it; loosen the residual tolerance given in !SOLVER'
        else
          msg(1) = '### WARNING:'//trim(noconv)
          msg(2) = '   the eigenvalues may carry an error of a comparable order'
        endif
        if(myrank == 0)then
          write(*,'(a/a)')    trim(msg(1)), trim(msg(2))
          write(ilog,'(a/a)') trim(msg(1)), trim(msg(2))
        endif
        if(unusable) call fstr_abort( hecmw_exit_noconv )
      endif
 
      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