fstr_solve_QuasiNewton.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_quasinewton
 
  use m_fstr_nonlinearmethod
  use m_fstr_iterationcontrol
  use m_fstr_nodalkinematics, only: fstr_apply_solution_increment, fstr_commit_solution_increment
 
  implicit none
  ! parameters for line search
  real(kind=kreal), parameter :: c_line_search=2.0, psi0_line_search=1.0d-2
  real(kind=kreal), parameter :: delta_wolfe=0.2, sigma_wolfe=0.9, eps_wolfe=1.0d-3
  real(kind=kreal), parameter :: omega_wolfe=0.001, delta_approx_wolfe=0.7
  real(kind=kreal) :: c_wolfe, q_wolfe
  integer, parameter :: n_mem_max=10
  integer(kind=kint), parameter :: maxiter_ls = 10
  integer(kind=kint), parameter :: pot_type=1
 
contains
 
  subroutine fstr_quasi_newton( cstep, hecMESH, hecMAT, fstrSOLID, fstrPARAM, &
    restrt_step_num, sub_step, ctime, dtime )
    implicit none
 
    integer, intent(in)                   :: cstep
    type (hecmwST_local_mesh)             :: hecMESH
    type (hecmwST_matrix)                 :: hecMAT
    type (fstr_solid)                     :: fstrSOLID
    integer, intent(in)                   :: sub_step
    real(kind=kreal), intent(in)          :: ctime     
    real(kind=kreal), intent(in)          :: dtime     
    type (fstr_param)                     :: fstrPARAM
    type (hecmwST_matrix_lagrange)        :: hecLagMAT
 
    type (hecmwST_local_mesh), pointer :: hecMESHmpc
    type (hecmwST_matrix), pointer :: hecMATmpc
    integer(kind=kint) :: ndof
    integer(kind=kint) :: i, iter
    integer(kind=kint) :: stepcnt
    integer(kind=kint) :: restrt_step_num
    real(kind=kreal)   :: tt0, tt, res, qnrm, rres, tincr, xnrm, dunrm, rxnrm
    logical :: isLinear = .false.
    integer(kind=kint) :: iterStatus
 
    real(kind=kreal), allocatable :: z_k(:), s_k(:,:), y_k(:,:), g_prev(:), rho_k(:)
    real(kind=kreal) :: sdoty
    integer :: n_mem
    integer :: len_vector
    integer(kind=kint) :: k
 
    integer :: u_debug
    integer :: max_iter_bak
    logical :: flag_approx_Wolfe
 
    max_iter_bak = fstrsolid%step_ctrl(cstep)%max_iter
    fstrsolid%step_ctrl(cstep)%max_iter = 100*fstrsolid%step_ctrl(cstep)%max_iter
 
    call hecmw_mpc_mat_init(hecmesh, hecmat, hecmeshmpc, hecmatmpc)
 
    if(.not. fstrpr%nlgeom)then
      islinear = .true.
    endif
 
    hecmat%NDOF = hecmesh%n_dof
    ndof = hecmat%NDOF
 
    stepcnt = 0
 
    tincr = dtime
    if( fstrsolid%step_ctrl(cstep)%solution == stepstatic ) tincr = 0.d0
    call fstr_init_newton(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, heclagmat, ndof, 0, hecmat)
    fstrsolid%GL0(:) = fstrsolid%GL(:) !store external load at du=0
 
    !! initialize du for non-zero Dirichlet condition
    call fstr_addbc(cstep, hecmesh, hecmat, fstrsolid, fstrparam, heclagmat, 1, rhsvector=fstrsolid%dunode)
    !! update residual vector
    call fstr_calc_residual_vector(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam)
 
    len_vector = hecmesh%n_node*ndof
    allocate(z_k(len_vector))
    allocate(s_k(len_vector, n_mem_max))
    allocate(y_k(len_vector, n_mem_max))
    allocate(g_prev(len_vector))
    allocate(rho_k(n_mem_max))
    z_k(:) = 0.0d0
    s_k(:,:) = 0.0d0
    y_k(:,:) = 0.0d0
    g_prev(:) = 0.0d0
    rho_k(:) = 0.0d0
    do i=1,len_vector
      y_k(i,1) = -hecmat%B(i)
    enddo
 
    ! parameter to judge Wolfe/approx Wolfe selection
    c_wolfe = 0.0d0
    q_wolfe = 0.0d0
    flag_approx_wolfe = .false.
    n_mem = 1
    ! ----- Inner Iteration, lagrange multiplier constant
    do iter=1,fstrsolid%step_ctrl(cstep)%max_iter
      stepcnt = stepcnt+1
 
      ! ----- calculate search direction by limited BFGS method
      do i=1,hecmesh%n_node*ndof
          g_prev(i) = -hecmat%B(i)
      enddo
      call fstr_calc_direction_lbfgs(hecmesh, g_prev, s_k, y_k, rho_k, z_k, n_mem)
 
      ! ! ----- Set Boundary condition
      call fstr_addbc_to_direction_vector(z_k, hecmesh,fstrsolid, cstep)
 
      !----- line search of step length
      call fstr_line_search_along_direction(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k)
 
      ! ----- update the small displacement and the displacement for 1step
      !       \delta u^k => solver's solution
      !       \Delta u_{n+1}^{k} = \Delta u_{n+1}^{k-1} + \delta u^k
      call fstr_apply_solution_increment( hecmesh, fstrsolid, ndof, hecmat%X )
 
      !! set du for non-zero Dirichlet condition
      ! call fstr_AddBC(cstep, hecMESH, hecMAT, fstrSOLID, fstrPARAM, hecLagMAT, 1, RHSvector=fstrSOLID%dunode)
      call fstr_calc_residual_vector(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam)
 
      ! ----- check convergence
      call fstr_check_convergence(hecmesh, hecmat, fstrsolid, fstrpr, &
          ndof, iter, sub_step, cstep, &
          hecmat%B, 0, &
          res, res, &
          0, &
          iterstatus)
      if (iterstatus == kitrconverged) exit
      if (iterstatus == kitrdiverged .or. iterstatus==kitrfloatingerror) return
      ! if (iterStatus == kitrDiverged) exit
      ! if (iterStatus==kitrFloatingError) return
 
      n_mem = min(n_mem+1, n_mem_max)
      do k=n_mem, 2, -1
        s_k(:,k) = s_k(:,k-1)
        y_k(:,k) = y_k(:,k-1)
        rho_k(k) = rho_k(k-1)
      enddo
 
      do i = 1, hecmesh%n_node*ndof
        s_k(i,1) = hecmat%X(i)
        y_k(i,1) = -hecmat%B(i) - g_prev(i)
      enddo
      call hecmw_innerproduct_r(hecmesh,ndof,s_k(:,1), y_k(:,1), sdoty)
      if (abs(sdoty) < 1.0d-10) then
        rho_k(1) = 0.0d0
      else
        rho_k(1) = 1.0d0/sdoty
      endif
    enddo
    ! ----- end of inner loop
 
    fstrsolid%NRstat_i(knstmaxit) = max(fstrsolid%NRstat_i(knstmaxit),iter) ! logging newton iteration(maxtier)
    fstrsolid%NRstat_i(knstsumit) = fstrsolid%NRstat_i(knstsumit) + iter    ! logging newton iteration(sum of iter)
 
    ! ----- update the total displacement
    ! u_{n+1} = u_{n} + \Delta u_{n+1}
      call fstr_commit_solution_increment( hecmesh, fstrsolid, ndof )
 
    call fstr_updatestate( hecmesh, fstrsolid, tincr )
 
    fstrsolid%CutBack_stat = 0
    call hecmw_mpc_mat_finalize(hecmesh, hecmat, hecmeshmpc, hecmatmpc)
 
    fstrsolid%step_ctrl(cstep)%max_iter = max_iter_bak
  end subroutine fstr_quasi_newton
 
  subroutine fstr_calc_direction_lbfgs(hecMESH, g_prev, s_k, y_k, rho_k, z_k, n_mem)
    implicit none
 
    type (hecmwST_local_mesh)             :: hecMESH
    real(kind=kreal) :: z_k(:), s_k(:,:), y_k(:,:), g_prev(:), rho_k(:)
    integer :: n_mem
 
    real(kind=kreal), allocatable :: q(:)
    real(kind=kreal) :: alpha(n_mem), beta
    real(kind=kreal) :: sdotq, ysq, gamma, ydotz, g_max
    
    integer :: len_vector, ndof
    integer(kind=kint) :: k,i
 
    ndof = hecmesh%n_dof
    len_vector = hecmesh%n_node*hecmesh%n_dof
    allocate(q(len_vector))
 
    q(1:len_vector) = g_prev(1:len_vector)
 
    do k=1, n_mem
      call hecmw_innerproduct_r(hecmesh,ndof,s_k(:,k), q, sdotq)
      alpha(k) = rho_k(k) * sdotq
      do i=1, len_vector
        q(i) = q(i) - alpha(k)*y_k(i,k)
      enddo
    enddo
    call hecmw_innerproduct_r(hecmesh,ndof,y_k(:,1), y_k(:,1), ysq)
    if (n_mem==1) then
      call hecmw_absmax_r(hecmesh, ndof, g_prev, g_max)
      ! if (g_max==0.0d0) then
      !   write(6,*) 'gradient of potential is zero-vector'
      !   stop
      ! endif
      ! gamma = 1.0d0/g_max
      gamma = 1.0d0
    else if (abs(rho_k(1)) < 1.0d-10) then
      gamma = 1.0d0
    else
      gamma = 1.0d0/(rho_k(1)*ysq)
    endif
 
    do i=1, len_vector
      z_k(i) = gamma*q(i)
    enddo
 
    do k=n_mem, 1, -1
      call hecmw_innerproduct_r(hecmesh,ndof,y_k(:,k), z_k, ydotz)
      beta = rho_k(k)*ydotz
      do i=1, len_vector
        z_k(i)=z_k(i) + s_k(i,k)*(alpha(k)-beta)
      enddo
    enddo
    deallocate(q)
  end subroutine fstr_calc_direction_lbfgs
 
  subroutine fstr_addbc_to_direction_vector(z_k, hecMESH,fstrSOLID, cstep)
    implicit none
    type (hecmwST_local_mesh)             :: hecMESH
    type (fstr_solid)                     :: fstrSOLID
    integer(kind=kint) :: cstep
    real(kind=kreal) :: z_k(:)
 
    integer(kind=kint) :: ig0, grpid, ig, ityp, idofS, idofE, iS0, iE0, ik, in, idof, ndof
 
    ndof = hecmesh%n_dof
    !   ----- Prescibed displacement Boundary Conditions
    do ig0 = 1, fstrsolid%BOUNDARY_ngrp_tot
      grpid = fstrsolid%BOUNDARY_ngrp_GRPID(ig0)
      if( .not. fstr_isboundaryactive( fstrsolid, grpid, cstep ) ) cycle
      ig   = fstrsolid%BOUNDARY_ngrp_ID(ig0)
      ityp = fstrsolid%BOUNDARY_ngrp_type(ig0)
      idofs = ityp/10
      idofe = ityp - idofs*10
      !
      is0 = hecmesh%node_group%grp_index(ig-1) + 1
      ie0 = hecmesh%node_group%grp_index(ig  )
      !
      do ik = is0, ie0
        in = hecmesh%node_group%grp_item(ik)
        !
        do idof = idofs, idofe
            z_k(ndof*(in-1)+idof) = 0.0d0
        enddo
      enddo
    enddo
  end subroutine fstr_addbc_to_direction_vector
 
  subroutine fstr_apply_alpha0(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM, z_k, h_prime, pot)
    implicit none
    type (hecmwST_local_mesh)             :: hecMESH
    type (hecmwST_matrix)                 :: hecMAT
    type (fstr_solid)                     :: fstrSOLID
    real(kind=kreal), intent(in)          :: ctime     
    real(kind=kreal), intent(in) :: tincr
    integer(kind=kint) :: iter
    integer, intent(in)                   :: cstep
    real(kind=kreal), intent(in)          :: dtime     
    type (fstr_param)                     :: fstrPARAM
    real(kind=kreal), intent(in) :: z_k(:)
    real(kind=kreal) :: h_prime, pot
 
    hecmat%X(:) = 0.0d0
    call hecmw_innerproduct_r(hecmesh, hecmat%NDOF, hecmat%B, z_k, h_prime)
    pot = fstr_get_potential(cstep,hecmesh,hecmat,fstrsolid,pot_type)
  end subroutine fstr_apply_alpha0
 
  subroutine fstr_apply_alpha(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM, z_k, alpha, h_prime, pot)
    implicit none
    type (hecmwST_local_mesh)             :: hecMESH
    type (hecmwST_matrix)                 :: hecMAT
    type (fstr_solid)                     :: fstrSOLID
    real(kind=kreal), intent(in)          :: ctime     
    real(kind=kreal), intent(in) :: tincr
    integer(kind=kint) :: iter
    integer, intent(in)                   :: cstep
    real(kind=kreal), intent(in)          :: dtime     
    type (fstr_param)                     :: fstrPARAM
    real(kind=kreal), intent(in) :: z_k(:)
    real(kind=kreal), intent(in) :: alpha
    real(kind=kreal) :: h_prime, pot
 
    hecmat%X(:) = -alpha*z_k(:)
    call fstr_calc_residual_vector_with_x(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam)
    call hecmw_innerproduct_r(hecmesh, hecmat%NDOF, hecmat%B, z_k, h_prime)
    pot = fstr_get_potential_with_x(cstep,hecmesh,hecmat,fstrsolid,pot_type)
  end subroutine fstr_apply_alpha
 
  subroutine fstr_init_line_search_range(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM, z_k, &
     h_prime_0, pot_0, alpha_S, h_prime_S, pot_S,  alpha_E, h_prime_E, pot_E)
    implicit none
    type (hecmwST_local_mesh)             :: hecMESH
    type (hecmwST_matrix)                 :: hecMAT
    type (fstr_solid)                     :: fstrSOLID
    real(kind=kreal), intent(in)          :: ctime     
    real(kind=kreal), intent(in) :: tincr
    integer(kind=kint) :: iter
    integer, intent(in)                   :: cstep
    real(kind=kreal), intent(in)          :: dtime     
    type (fstr_param)                     :: fstrPARAM
    real(kind=kreal), intent(in) :: z_k(:)
    real(kind=kreal), intent(in) :: h_prime_0, pot_0
    real(kind=kreal) :: alpha_s, h_prime_s, pot_s
    real(kind=kreal) :: alpha_e, h_prime_e, pot_e
 
    real(kind=kreal) :: alpha_s_new, h_prime_s_new, pot_s_new
    real(kind=kreal) :: alpha_e_new, h_prime_e_new, pot_e_new
    real(kind=kreal) :: alpha_tmp, h_prime_tmp, pot_tmp
    real(kind=kreal) :: z_max
    real(kind=kreal) :: pot_eps
    pot_eps = eps_wolfe*c_wolfe
 
    alpha_s = 0.0d0
    h_prime_s = h_prime_0
    pot_s = pot_0
 
    if (iter==1) then
      call hecmw_absmax_r(hecmesh, hecmat%NDOF, z_k, z_max)
      alpha_e = psi0_line_search/z_max
    else
      alpha_e = 1.0d0
    end if
    call fstr_apply_alpha(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, alpha_e, h_prime_e, pot_e)
 
    do while (h_prime_e < 0.0d0)
      ! if (pot_E <= pot_0 + pot_eps) then
        ! alpha_S = alpha_E
        ! h_prime_S = h_prime_E ! so h_prime_S < 0
        ! pot_S = h_prime_S
 
        alpha_e = alpha_e * c_line_search
        call fstr_apply_alpha(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, &
                                     &  fstrparam, z_k, alpha_e, h_prime_e, pot_e)
      ! else
      !   alpha_tmp = 2.0d0*alpha_E
      !   h_prime_tmp = 0.0d0 ! elemact value
      !   pot_tmp = 0.0d0 ! elemact value
      !   call fstr_get_new_range_with_potential(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM, z_k, pot_0, &
      !     alpha_S, h_prime_S, pot_S, alpha_tmp, h_prime_tmp, pot_tmp, alpha_E, h_prime_E, pot_E, &
      !     alpha_S_new, h_prime_S_new, pot_S_new, alpha_E_new, h_prime_E_new, pot_E_new)
 
      !   alpha_S = alpha_S_new
      !   h_prime_S = h_prime_S_new
      !   pot_S = pot_S_new
 
      !   alpha_E = alpha_E_new
      !   h_prime_E = h_prime_E_new
      !   pot_E = pot_E_new
      !   return
      ! end if
    enddo
  end subroutine fstr_init_line_search_range
 
  subroutine fstr_update_range(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM, z_k, &
     alpha_S, h_prime_S, pot_S,  alpha_E, h_prime_E, pot_E,  alpha_c, h_prime_c, pot_c, h_prime_0, pot_0)
    implicit none
    type (hecmwST_local_mesh)             :: hecMESH
    type (hecmwST_matrix)                 :: hecMAT
    type (fstr_solid)                     :: fstrSOLID
    real(kind=kreal), intent(in)          :: ctime     
    real(kind=kreal), intent(in) :: tincr
    integer(kind=kint) :: iter
    integer, intent(in)                   :: cstep
    real(kind=kreal), intent(in)          :: dtime     
    type (fstr_param)                     :: fstrPARAM
    real(kind=kreal) :: z_k(:)
    real(kind=kreal) :: alpha_s, h_prime_s, pot_s
    real(kind=kreal) :: alpha_e, h_prime_e, pot_e
    real(kind=kreal) :: alpha_c, h_prime_c, pot_c
    real(kind=kreal) :: h_prime_0, pot_0
 
    real(kind=kreal) :: alpha_a, h_prime_a, pot_a
    real(kind=kreal) :: alpha_b, h_prime_b, pot_b
    real(kind=kreal) :: alpha_a_bar, h_prime_a_bar, pot_a_bar
    real(kind=kreal) :: alpha_b_bar, h_prime_b_bar, pot_b_bar
    real(kind=kreal) :: alpha_c_bar, h_prime_c_bar, pot_c_bar
 
    call fstr_get_new_range_with_potential(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, pot_0, &
      alpha_s, h_prime_s, pot_s, alpha_e, h_prime_e, pot_e, alpha_c, h_prime_c, pot_c, &
      alpha_a, h_prime_a, pot_a, alpha_b, h_prime_b, pot_b)
 
    if (alpha_c == alpha_a) then
      call fstr_get_secant(alpha_s, h_prime_s, alpha_a, h_prime_a, alpha_c_bar)
      call fstr_apply_alpha(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, &
        alpha_c_bar, h_prime_c_bar, pot_c_bar)
    end if
    if (alpha_c == alpha_b) then
      call fstr_get_secant(alpha_b, h_prime_b, alpha_e, h_prime_e, alpha_c_bar)
      call fstr_apply_alpha(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, &
        alpha_c_bar, h_prime_c_bar, pot_c_bar)
    end if
 
    if ((alpha_c == alpha_a) .or. (alpha_c == alpha_b)) then
      call fstr_get_new_range_with_potential(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, pot_0, &
        alpha_a, h_prime_a, pot_a, alpha_b, h_prime_b, pot_b, alpha_c_bar, h_prime_c_bar, pot_c_bar, &
        alpha_a_bar, h_prime_a_bar, pot_a_bar, alpha_b_bar, h_prime_b_bar, pot_b_bar)
    else
      alpha_a_bar = alpha_a
      h_prime_a_bar = h_prime_a
      pot_a_bar = pot_a
 
      alpha_b_bar = alpha_b
      h_prime_b_bar = h_prime_b
      pot_b_bar = pot_b
    end if
 
 
    if ((alpha_b_bar - alpha_a_bar) > 0.66*(alpha_e - alpha_s)) then
      alpha_c_bar = 0.5*(alpha_a_bar+alpha_b_bar)
      call fstr_apply_alpha(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, &
        alpha_c_bar, h_prime_c_bar, pot_c_bar)
      call fstr_get_new_range_with_potential(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, pot_0, &
        alpha_a, h_prime_a, pot_a, alpha_b, h_prime_b, pot_b, alpha_c_bar, h_prime_c_bar, pot_c_bar, &
        alpha_s, h_prime_s, pot_s, alpha_e, h_prime_e, pot_e)
    else
      alpha_s = alpha_a_bar
      h_prime_s = h_prime_a_bar
      pot_s = pot_a_bar
 
      alpha_e = alpha_b_bar
      h_prime_e = h_prime_b_bar
      pot_e = pot_b_bar
    end if
  end subroutine fstr_update_range
 
  subroutine fstr_get_new_range_with_potential(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, &
    fstrPARAM, z_k, pot_0, alpha_a, h_prime_a, pot_a, alpha_b, h_prime_b, pot_b, alpha_c, h_prime_c, pot_c, &
     alpha_a_bar, h_prime_a_bar, pot_a_bar, alpha_b_bar, h_prime_b_bar, pot_b_bar)
    implicit none
    type (hecmwST_local_mesh)             :: hecMESH
    type (hecmwST_matrix)                 :: hecMAT
    type (fstr_solid)                     :: fstrSOLID
    real(kind=kreal), intent(in)          :: ctime     
    real(kind=kreal), intent(in) :: tincr
    integer(kind=kint) :: iter
    integer, intent(in)                   :: cstep
    real(kind=kreal), intent(in)          :: dtime     
    type (fstr_param)                     :: fstrPARAM
    real(kind=kreal), intent(in) :: z_k(:)
    real(kind=kreal), intent(in) :: pot_0
    real(kind=kreal), intent(in) :: alpha_a, h_prime_a, pot_a
    real(kind=kreal), intent(in) :: alpha_b, h_prime_b, pot_b
    real(kind=kreal), intent(in) :: alpha_c, h_prime_c, pot_c
    real(kind=kreal) :: alpha_a_bar, h_prime_a_bar, pot_a_bar
    real(kind=kreal) :: alpha_b_bar, h_prime_b_bar, pot_b_bar
 
    integer, parameter :: count_max=100
    integer :: count_while
    real(kind=kreal) :: alpha_d, h_prime_d, pot_d
    real(kind=kreal) :: theta_ls = 0.5d0
    real(kind=kreal) :: pot_eps
    pot_eps = eps_wolfe*c_wolfe
 
    ! case of NOT (a < c < b)
    if (alpha_c <= alpha_a .or. alpha_b <= alpha_c) then
      alpha_a_bar = alpha_a
      h_prime_a_bar = h_prime_a
      pot_a_bar = pot_a
 
      alpha_b_bar = alpha_b
      h_prime_b_bar = h_prime_b
      pot_b_bar = pot_b
      return
    end if
 
    if (h_prime_c >= 0.0d0) then
      alpha_a_bar = alpha_a
      h_prime_a_bar = h_prime_a
      pot_a_bar = pot_a
 
      alpha_b_bar = alpha_c
      h_prime_b_bar = h_prime_c
      pot_b_bar = pot_c
      return
    end if
    ! below here, it can be assumed that h_prime_c<0
 
    if(pot_c <= pot_0 + pot_eps) then
      alpha_a_bar = alpha_c
      h_prime_a_bar = h_prime_c
      pot_a_bar = pot_c
 
      alpha_b_bar = alpha_b
      h_prime_b_bar = h_prime_b
      pot_b_bar = pot_b
      return
    end if
 
    alpha_a_bar = alpha_a
    h_prime_a_bar = h_prime_a
    pot_a_bar = pot_a
 
    alpha_b_bar = alpha_b
    h_prime_b_bar = h_prime_b
    pot_b_bar = pot_b
 
    count_while = 0
    do while(count_while < count_max)
      count_while = count_while + 1
 
      alpha_d = (1.0d0-theta_ls)*alpha_a_bar + theta_ls*alpha_b_bar
      call fstr_apply_alpha(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, alpha_d, h_prime_d, pot_d)
 
      if (h_prime_d >= 0.0d0) then
        alpha_b_bar = alpha_d
        h_prime_b_bar = h_prime_d
        pot_b_bar = pot_d
        return
      end if
 
      if(pot_d <= pot_0 + pot_eps) then
        alpha_a_bar = alpha_d
        h_prime_a_bar = h_prime_d
        pot_a_bar = pot_d
      else
        alpha_b_bar = alpha_d
        h_prime_b_bar = h_prime_d
        pot_b_bar = pot_d
      end if
    end do
    write(6,*) 'fstr_get_new_range_with_potential reached loop count max.', hecmesh%my_rank, alpha_a_bar, alpha_b_bar
  end subroutine fstr_get_new_range_with_potential
 
 
  subroutine fstr_line_search_along_direction(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM, z_k)
    implicit none
    type (hecmwST_local_mesh)             :: hecMESH
    type (hecmwST_matrix)                 :: hecMAT
    type (fstr_solid)                     :: fstrSOLID
    real(kind=kreal), intent(in)          :: ctime     
    real(kind=kreal), intent(in) :: tincr
    integer(kind=kint) :: iter
    integer, intent(in)                   :: cstep
    real(kind=kreal), intent(in)          :: dtime     
    type (fstr_param)                     :: fstrPARAM
    real(kind=kreal) :: z_k(:)
 
    real(kind=kreal) :: alpha_s, h_prime_s, pot_s
    real(kind=kreal) :: alpha_e, h_prime_e, pot_e
    real(kind=kreal) :: alpha_c, h_prime_c, pot_c
    real(kind=kreal) :: h_prime_0, pot_0
    logical :: flag_converged
    integer :: ndof, len_vector
    real(kind=kreal) :: res
 
    integer(kind=kint) :: i, ierr, iter_ls
    real(kind=kreal) :: z_max
    integer :: elemact
 
    ndof = hecmat%NDOF
    len_vector = hecmesh%n_node*hecmesh%n_dof
 
    call fstr_apply_alpha0(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, h_prime_0, pot_0)
    if (h_prime_0 > 0.0d0) then
      write(6,*) 'residual vector is not directed to potential decretion.', h_prime_0
      stop
    endif
 
    call fstr_init_line_search_range(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, &
    fstrparam, z_k, h_prime_0, pot_0, alpha_s, h_prime_s, pot_s,  alpha_e, h_prime_e, pot_e)
    if( hecmesh%my_rank == 0 ) then
      write(6,'(a, 6es27.16e3)') 'range initialized: alpha_S, alpha_E, h_prime_S, h_prime_E, pot_S, pot_E', &
          &  alpha_s, alpha_e, h_prime_s, h_prime_e, pot_s, pot_e
    endif    
 
        q_wolfe = 1 + delta_approx_wolfe * q_wolfe
        c_wolfe = c_wolfe + (abs(pot_0)-c_wolfe) / q_wolfe
 
    flag_converged = .false.
    iter_ls=0
    do while (iter_ls<maxiter_ls)
      call fstr_get_secant(alpha_s, h_prime_s, alpha_e, h_prime_e, alpha_c)
      call fstr_apply_alpha(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, alpha_c, h_prime_c, pot_c)
 
      if (abs( pot_c - pot_0) <= (omega_wolfe * c_wolfe) ) then
        flag_converged = fstr_approx_wolfe_condition(alpha_c, h_prime_c, pot_c, h_prime_0, pot_0)
      else
        flag_converged = fstr_wolfe_condition(alpha_c, h_prime_c, pot_c, h_prime_0, pot_0)
      endif
      if (flag_converged) exit
 
      call fstr_update_range(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, &
        alpha_s, h_prime_s, pot_s,  alpha_e, h_prime_e, pot_e,  alpha_c, h_prime_c, pot_c, h_prime_0, pot_0)
      iter_ls = iter_ls +1
    enddo
    if( hecmesh%my_rank == 0 ) then
      write(6,'(a, 1i8, 5es27.16e3)') 'converged: alpha_S, alpha_E, alpha_c, h_prime_c, pot_c', &
          &  iter_ls, &
          &  alpha_s, alpha_e, alpha_c, h_prime_c, pot_c
    endif    
  end subroutine fstr_line_search_along_direction
 
  subroutine fstr_get_secant(a, Fa, b, Fb, c)
    implicit none
    real(kind=kreal), intent(in) :: a,b, fa,fb
    real(kind=kreal), intent(out) :: c
    if(fb /= fa) then
      c = (a*fb  - b*fa) / ( fb - fa )
    else
      c = 0.5*a + 0.5*b
    endif
  end subroutine fstr_get_secant
 
  function fstr_wolfe_condition(alpha_c, h_prime_c, pot_c, h_prime_0, pot_0) result(flag_converged)
    implicit none
    logical :: flag_converged
 
    real(kind=kreal), intent(in) :: alpha_c, h_prime_0, h_prime_c, pot_0, pot_c
    real(kind=kreal) :: wolfe1_left, wolfe1_right, wolfe2_left, wolfe2_right
 
    wolfe1_left = pot_c - pot_0
    wolfe1_right = delta_wolfe*h_prime_0*alpha_c
 
    wolfe2_left = h_prime_c
    wolfe2_right = sigma_wolfe*h_prime_0
 
    flag_converged = (wolfe1_left<=wolfe1_right) .and. (wolfe2_left >= wolfe2_right)
    if (flag_converged) write(6,'(a, 4es27.16e3)') 'oWolfe: ', wolfe1_left, wolfe1_right, wolfe2_left, wolfe2_right
  end function fstr_wolfe_condition
 
  function fstr_approx_wolfe_condition(alpha_c, h_prime_c, pot_c, h_prime_0, pot_0) result(flag_converged)
    implicit none
    logical :: flag_converged
 
    real(kind=kreal), intent(in) :: alpha_c, h_prime_0, h_prime_c, pot_0, pot_c
    real(kind=kreal) :: wolfe1_left, wolfe1_right, wolfe2_left, wolfe2_right
 
    real(kind=kreal) :: pot_eps
    pot_eps = eps_wolfe*c_wolfe
 
    wolfe1_left =  ( 2.0 * delta_wolfe - 1.0d0 ) * h_prime_0
    wolfe1_right = h_prime_c
 
    wolfe2_left = h_prime_c
    wolfe2_right = sigma_wolfe*h_prime_0
 
    flag_converged = &
      (wolfe1_left>=wolfe1_right) &
      .and. (wolfe2_left >= wolfe2_right) &
      .and. (pot_c <= pot_0 + pot_eps)
    if (flag_converged) write(6,'(a, 6es27.16e3)') 'aWolfe: ', wolfe1_left, wolfe1_right, &
    &  wolfe2_left, wolfe2_right, pot_c, pot_0 + (eps_wolfe*c_wolfe)
  end function fstr_approx_wolfe_condition
end module m_fstr_quasinewton