text:
module unstructured_grid_mod
implicit none
private
public :: icon_icosohedral_grid_constructor

type, public, abstract :: unstructured_grid
  private
  integer, public :: num_points
  integer, dimension(:,:), pointer :: cell_neighbors
  integer, dimension(:,:), pointer :: cell_secondary_neighbors
  integer :: num_primary_neighbors
  integer :: num_secondary_neighbors
  contains
    private
      procedure         :: initialise_unstructured_grid
      procedure, public :: get_cell_secondary_neighbors
      procedure, public :: generate_cell_neighbors
      procedure, public :: generate_cell_secondary_neighbors
      procedure, public :: generate_edge_cells
      procedure, public :: generate_subfield_indices
      procedure, public :: generate_full_field_indices
      procedure, public :: unstructured_grid_destructor
      procedure(calculate_secondary_neighbors), deferred, public :: &
        calculate_secondary_neighbors
end type unstructured_grid

abstract interface
  subroutine calculate_secondary_neighbors(this)
    import unstructured_grid
    class(unstructured_grid), intent(inout) :: this
  end subroutine calculate_secondary_neighbors
end interface

type, public, extends(unstructured_grid) :: icon_icosohedral_grid
  private
  contains
    private
      procedure :: initialise_icon_icosohedral_grid
      procedure, public :: calculate_secondary_neighbors => icon_icosohedral_calculate_secondary_neighbors
end type icon_icosohedral_grid

interface icon_icosohedral_grid
  procedure :: icon_icosohedral_grid_constructor
end interface icon_icosohedral_grid

contains

  function get_cell_secondary_neighbors(this) result(secondary_neighbors)
    class(unstructured_grid), intent(inout) :: this
    integer, dimension(:,:), pointer :: secondary_neighbors
      secondary_neighbors => this%cell_secondary_neighbors
  end function get_cell_secondary_neighbors

  subroutine initialise_unstructured_grid(this,cell_neighbors_in, &
                                          num_primary_neighbors_in, &
                                          num_secondary_neighbors_in)
    class(unstructured_grid), intent(inout) :: this
    integer, dimension(:,:), pointer, intent(in) :: cell_neighbors_in
    integer, intent(in) :: num_primary_neighbors_in
    integer, intent(in) :: num_secondary_neighbors_in
      this%num_points = size(cell_neighbors_in,1)
      this%cell_neighbors => cell_neighbors_in
      this%num_primary_neighbors = num_primary_neighbors_in
      this%num_secondary_neighbors = num_secondary_neighbors_in
  end subroutine initialise_unstructured_grid

  subroutine unstructured_grid_destructor(this)
    class(unstructured_grid), intent(inout) :: this
      deallocate(this%cell_secondary_neighbors)
  end subroutine unstructured_grid_destructor

  subroutine initialise_icon_icosohedral_grid(this,cell_neighbors_in)
    class(icon_icosohedral_grid), intent(inout) :: this
    integer, dimension(:,:), pointer, intent(in) :: cell_neighbors_in
    integer :: num_primary_neighbors
    integer :: num_secondary_neighbors
      num_primary_neighbors = 3
      num_secondary_neighbors = 9
      call this%initialise_unstructured_grid(cell_neighbors_in, &
                                             num_primary_neighbors, &
                                             num_secondary_neighbors)
  end subroutine initialise_icon_icosohedral_grid

  function icon_icosohedral_grid_constructor(cell_neighbors_in) result(constructor)
    integer, dimension(:,:), pointer, intent(in) :: cell_neighbors_in
    type(icon_icosohedral_grid), pointer :: constructor
      allocate(constructor)
      call constructor%initialise_icon_icosohedral_grid(cell_neighbors_in)
  end function icon_icosohedral_grid_constructor

  function generate_edge_cells(this,cell_neighbors,cell_secondary_neighbors) &
      result(edge_cells)
    class(unstructured_grid), intent(in) :: this
    integer, dimension(:,:), pointer, intent(in) :: cell_neighbors
    integer, dimension(:,:), pointer, intent(in) :: cell_secondary_neighbors
    integer, dimension(:), pointer :: edge_cells
    integer, dimension(:), pointer :: expanded_edge_cells
    integer :: num_points_subarray
    integer :: num_edge_cells
    integer :: i,j
    logical :: is_edge_cell
      num_points_subarray = size(cell_neighbors,1)
      allocate(expanded_edge_cells(num_points_subarray))
      num_edge_cells = 0
      do i = 1,num_points_subarray
        is_edge_cell = .false.
        do j = 1,this%num_primary_neighbors
          if(cell_neighbors(i,j) == -1) then
            num_edge_cells = num_edge_cells + 1
            expanded_edge_cells(num_edge_cells) = i
            is_edge_cell = .true.
            exit
          end if
        end do
        if (is_edge_cell) exit
        do j = 1,this%num_secondary_neighbors
          if(cell_secondary_neighbors(i,j) == -1) then
            num_edge_cells = num_edge_cells + 1
            expanded_edge_cells(num_edge_cells) = i
            is_edge_cell = .true.
            exit
          end if
        end do
        if (is_edge_cell) exit
      end do
      allocate(edge_cells(num_edge_cells))
      do i = 1,num_edge_cells
        edge_cells(i) = expanded_edge_cells(i)
      end do
  end function generate_edge_cells

  function generate_full_field_indices(this,mask,subfield_indices) result(full_field_indices)
      class(unstructured_grid), intent(in) :: this
      integer, dimension(:), pointer, intent(in) :: subfield_indices
      logical, dimension(:), pointer, intent(in) :: mask
      integer, dimension(:), pointer :: full_field_indices
      integer :: num_points_subarray
      integer :: i
        num_points_subarray = count(mask)
        allocate(full_field_indices(num_points_subarray))
        do i = 1,this%num_points
          if(mask(i)) then
            full_field_indices(subfield_indices(i)) = i
          end if
        end do
  end function generate_full_field_indices

  function generate_subfield_indices(this,mask) result(subfield_indices)
      class(unstructured_grid), intent(in) :: this
      logical, dimension(:), pointer, intent(in) :: mask
      integer, dimension(:), pointer :: subfield_indices
      integer :: i, i_new
        allocate(subfield_indices(this%num_points))
        i_new = 1
        do i = 1,this%num_points
          if(mask(i)) then
            subfield_indices(i) = i_new
            i_new = i + i_new
          else
            subfield_indices(i) = -1
          end if
        end do
  end function generate_subfield_indices

  function generate_cell_neighbors(this,mask,subfield_indices,full_field_indices) &
      result(cell_neighbors)
    class(unstructured_grid), intent(in) :: this
    logical, dimension(:), pointer, intent(in) :: mask
    integer, dimension(:), pointer, intent(in) :: subfield_indices
    integer, dimension(:), pointer, intent(in) :: full_field_indices
    integer, dimension(:,:), pointer :: cell_neighbors
    integer :: num_points_subarray
    integer :: i,nbr_num
      if(all(mask)) then
        cell_neighbors => this%cell_neighbors
      else
        num_points_subarray = count(mask)
        allocate(cell_neighbors(num_points_subarray,this%num_primary_neighbors))
        do i = 1,num_points_subarray
          do nbr_num = 1,this%num_primary_neighbors
            cell_neighbors(i,nbr_num) = subfield_indices(this%cell_neighbors(full_field_indices(i),nbr_num))
          end do
        end do
      end if
  end function generate_cell_neighbors

  function generate_cell_secondary_neighbors(this,subfield_indices,full_field_indices,cell_neighbors) &
      result(cell_secondary_neighbors)
    class(unstructured_grid), intent(inout) :: this
    integer, dimension(:,:), pointer, intent(in) :: cell_neighbors
    integer, dimension(:,:), pointer :: cell_secondary_neighbors
    integer, dimension(:), pointer, intent(in) :: subfield_indices
    integer, dimension(:), pointer, intent(in) :: full_field_indices
    integer :: num_points_subarray
    integer :: i,j
      num_points_subarray = size(cell_neighbors)
      if( .not. associated(this%cell_secondary_neighbors)) then
        call this%calculate_secondary_neighbors()
      end if
      allocate(cell_secondary_neighbors(num_points_subarray,this%num_secondary_neighbors))
      do i = 1,num_points_subarray
        do j = 1,this%num_secondary_neighbors
          cell_secondary_neighbors(i,j) = subfield_indices(this%cell_secondary_neighbors(full_field_indices(i),j))
        end do
      end do
  end function generate_cell_secondary_neighbors

  subroutine icon_icosohedral_calculate_secondary_neighbors(this)
    class(icon_icosohedral_grid), intent(inout) :: this
    integer :: index_over_grid
    integer :: index_over_primary_nbrs,index_over_secondary_nbrs
    integer :: index_over_tertiary_nbrs
    integer :: primary_neighbor_index,secondary_neighbor_index
    integer :: tertiary_neighbor_index
    integer :: second_tertiary_neighbor_index
    integer :: valid_secondary_nbr_count
    integer :: gap_index = 2
    integer :: first_secondary_neighbor_index,second_secondary_neighbor_index
    integer :: second_index_over_secondary_nbrs
    integer :: second_index_over_tertiary_nbrs
    integer :: no_neighbor
      no_neighbor = 0
      allocate(this%cell_secondary_neighbors(this%num_points,9))
      do index_over_grid=1,this%num_points
        !Six secondary neighbors are neighbors of primary neighbors
        do index_over_primary_nbrs=1,3
          primary_neighbor_index = this%cell_neighbors(index_over_grid,index_over_primary_nbrs)
          valid_secondary_nbr_count = 1;
          do index_over_secondary_nbrs=1,3
            !2 rather than 3 times primary neighbor index as we miss out 1 secondary neighbor for each
            !primary neighbor
            secondary_neighbor_index = this%cell_neighbors(primary_neighbor_index,index_over_secondary_nbrs)
            if (secondary_neighbor_index /= index_over_grid) then
              !Note this leaves gaps for the remaining three secondary neighbors
              this%cell_secondary_neighbors(index_over_grid,3*(index_over_primary_nbrs - 1) + &
                                            valid_secondary_nbr_count) = secondary_neighbor_index
              valid_secondary_nbr_count = valid_secondary_nbr_count + 1
            end if
          end do
        end do
        !Three secondary neighbors are common neighbors of the existing secondary neighbors
        gap_index = 3;
        !Last secondary neighbor is as yet unfilled so loop only up to an index of 8
        do index_over_secondary_nbrs=1,8
          !skip as yet unfilled entries in the secondary neighbors array
          if (mod(index_over_secondary_nbrs,3) == 0) cycle
          first_secondary_neighbor_index = &
              this%cell_secondary_neighbors(index_over_grid,index_over_secondary_nbrs)
          !Last secondary neighbor is as yet unfilled so loop only up to an index of 7
          do second_index_over_secondary_nbrs=index_over_secondary_nbrs+2,8
            if (mod(second_index_over_secondary_nbrs,3) == 0) cycle
              second_secondary_neighbor_index = &
                this%cell_secondary_neighbors(index_over_grid,second_index_over_secondary_nbrs)
            !Some tertiary neighbors are also secondary neighbors
            do index_over_tertiary_nbrs=1,3
              tertiary_neighbor_index = &
                this%cell_neighbors(first_secondary_neighbor_index,index_over_tertiary_nbrs)
              !Test to see if this one of the twelve 5-point vertices in the grid
              if (second_secondary_neighbor_index == tertiary_neighbor_index) then
                this%cell_secondary_neighbors(index_over_grid,gap_index) = no_neighbor
                gap_index = gap_index + 3;
                cycle;
              end if
              do second_index_over_tertiary_nbrs=1,3
                second_tertiary_neighbor_index = &
                  this%cell_neighbors(second_secondary_neighbor_index,second_index_over_tertiary_nbrs)
                if(second_tertiary_neighbor_index == tertiary_neighbor_index) then
                  this%cell_secondary_neighbors(index_over_grid,gap_index) = tertiary_neighbor_index
                  gap_index = gap_index + 3;
                end if
              end do
            end do
          end do
        end do
      end do
  end subroutine icon_icosohedral_calculate_secondary_neighbors

end module unstructured_grid_mod

max_stars_repo_name: ThomasRiddick/DynamicHD

max_issues_repo_name: ThomasRiddick/DynamicHD

max_forks_repo_name: ThomasRiddick/DynamicHD


!=====================================================================*
!                                                                     *
!   Software Name : ppohBEM                                           *
!         Version : 0.1                                               *
!                                                                     *
!   License                                                           *
!     This file is part of ppohBEM.                                   *
!     ppohBEM is a free software, you can use it under the terms   *
!     of The MIT License (MIT). See LICENSE file and User's guide     *
!     for more details.                                               *
!                                                                     *
!   ppOpen-HPC project:                                               *
!     Open Source Infrastructure for Development and Execution of     *
!     Large-Scale Scientific Applications on Post-Peta-Scale          *
!     Supercomputers with Automatic Tuning (AT).                      *
!                                                                     *
!   Organizations:                                                    *
!     The University of Tokyo                                         *
!       - Information Technology Center                               *
!       - Atmosphere and Ocean Research Institute (AORI)              *
!       - Interfaculty Initiative in Information Studies              *
!         /Earthquake Research Institute (ERI)                        *
!       - Graduate School of Frontier Science                         *
!     Kyoto University                                                *
!       - Academic Center for Computing and Media Studies             *
!     Hokkaido University                                             *
!       - Information Initiative Center                               *
!     Japan Agency for Marine-Earth Science and Technology (JAMSTEC)  *
!                                                                     *
!   Sponsorship:                                                      *
!     Japan Science and Technology Agency (JST), Basic Research       *
!     Programs: CREST, Development of System Software Technologies    *
!     for post-Peta Scale High Performance Computing.                 *
!                                                                     *
!   Copyright (c) 2012 <Takeshi Iwashita, Takeshi Mifune, Yuki Noseda,*
!                    Yasuhito Takahashi, Masatoshi Kawai, Akihiro Ida>*
!                                                                     *
!=====================================================================*
!
module m_ppohBEM_user_func
  integer, parameter :: KIND_REAL = kind(1d0)
  real(kind = KIND_REAL), parameter :: PI = 3.141592653589793
  real(kind = KIND_REAL), parameter :: EPSILON_0 = 8.854187818 * 1d-12

  public face_integral

  contains

    real(kind = KIND_REAL) function face_integral(xs, ys, zs, x, y, z)
    real(kind = KIND_REAL), intent(in) :: xs(3), ys(3), zs(3)
    real(kind = KIND_REAL), intent(in) :: x, y, z

      real(kind = KIND_REAL) :: r(3)
      real(kind = KIND_REAL) :: xi, xj, yi, dx, dy, t, l, m, d, ti, tj
      real(kind = KIND_REAL) :: theta, omega, q, g, zp, zpabs

      integer :: i, j
      real(kind = KIND_REAL) :: u(3), v(3), w(3)
      real(kind = KIND_REAL) :: ox, oy, oz

      r(:) = sqrt( (xs(:) - x)**2 + (ys(:) - y)**2 + (zs(:) - z)**2 )

      u(1) = xs(2) - xs(1);  v(1) = xs(3) - xs(2)
      u(2) = ys(2) - ys(1);  v(2) = ys(3) - ys(2)
      u(3) = zs(2) - zs(1);  v(3) = zs(3) - zs(2)
      call cross_product(u, v, w)
      w(:) = w(:) / sqrt( dot_product(w, w) )

      u(1) = x - xs(1);  u(2) = y - ys(1);  u(3) = z - zs(1)
      zp = dot_product(u, w)
      ox = x - zp * w(1);  oy = y - zp * w(2);  oz = z - zp * w(3)
      zpabs = abs(zp)

      face_integral = 0d0
      do i = 1, 3
         j = mod(i, 3) + 1
         u(1) = xs(j) - ox;  u(2) = ys(j) - oy; u(3) = zs(j) - oz;
               xj = sqrt( dot_product(u, u) )
         u(:) = u(:) / xj
         call cross_product(w, u, v)
         xi = (xs(i) - ox) * u(1) + (ys(i) - oy) * u(2) + (zs(i) - oz) * u(3)
         yi = (xs(i) - ox) * v(1) + (ys(i) - oy) * v(2) + (zs(i) - oz) * v(3)

         dx = xj - xi;  dy = - yi ! (yj = 0)
         t = sqrt( dx**2 + dy**2 )
         l = dx/t;  m = dy/t
         d = l * yi - m * xi
         ti = l * xi + m * yi;  tj = l * xj ! (yj = 0)

         !            theta = sign(1d0, yi) * acos( xi / sqrt( xi**2 + yi**2 ) ) ! (xj > 0, yj = 0)
         theta = atan2(yi, xi)
         omega = theta - atan2( r(i) * d, zpabs * ti ) + atan2( r(j) * d, zpabs * tj )
         q = log( (r(j) + tj) / ( r(i) + ti ) )
         g = d * q - zpabs * omega
         face_integral = face_integral + g
      enddo

      face_integral = abs(face_integral) / (4d0 * PI * EPSILON_0)

    end function

!***********************************************************************
    subroutine cross_product(u, v, w)
      real(kind = KIND_REAL) :: u(3), v(3), w(3)

      w(1) = u(2) * v(3) - u(3) * v(2)
      w(2) = u(3) * v(1) - u(1) * v(3)
      w(3) = u(1) * v(2) - u(2) * v(1)

    end subroutine cross_product

  end module m_ppohBEM_user_func


text:
module mod_convert_ubt_to_wbv
  use mod_prec
  use mod_error_id
  use mod_rotation
  use mod_orth_band_types
  use mod_band_types
  implicit none

  private

  public :: convert_ubt_to_wbv, d_convert_ubt_to_wbv, z_convert_ubt_to_wbv, &
       f_convert_ubt_to_wbv, f_d_convert_ubt_to_wbv, f_z_convert_ubt_to_wbv, &
       d_wbv_of_ubt, z_wbv_of_ubt, wbv_of

  interface convert_ubt_to_wbv
     module procedure d_convert_ubt_to_wbv, z_convert_ubt_to_wbv
  end interface convert_ubt_to_wbv

  interface f_convert_ubt_to_wbv
     module procedure f_d_convert_ubt_to_wbv, f_z_convert_ubt_to_wbv
  end interface f_convert_ubt_to_wbv

  interface wbv_of
     module procedure d_wbv_of_ubt, z_wbv_of_ubt
  end interface wbv_of

contains

  function d_wbv_of_ubt(ubt,error) result(wbv)
    type(d_wbv) :: wbv
    type(d_ubt), intent(in) :: ubt
    type(error_info), intent(inout), optional :: error

    type(d_ubt), allocatable :: ubt1
    type(routine_info), parameter :: info=info_d_wbv_of_ubt

    if (failure(error)) return
    call push_id(info, error)

    wbv=d_new_wbv(get_n(ubt), ubt%lbw, ubt%ubw)
    ubt1=d_new_ubt(get_n(ubt),ubt%lbw+1,  ubt%ubw+1)
    call copy(ubt1,ubt)
    call d_convert_ubt_to_wbv(ubt1,wbv,error)

    call pop_id(error)
  end function d_wbv_of_ubt

  ! Errors:
  ! 0: no error
  ! 1: n<1
  ! 2: Insufficient storage in ubt
  ! 3: Insufficient stroage in wbv
  ! 4: ubt%n /= wbv%n

  subroutine d_convert_ubt_to_wbv(ubt, wbv, error)
    type(d_ubt) :: ubt
    type(d_wbv) :: wbv
    type(error_info), intent(inout), optional :: error
    type(routine_info), parameter :: info=info_d_convert_ubt_to_wbv

    if (failure(error)) return
    call push_id(info, error)

    if (get_n(ubt) < 1) then
       call set_error(1, info, error); return
    end if
    if ((get_ubwmax(ubt) < ubt%ubw+1 .and. ubt%ubw < get_n(ubt)-1) .or. &
         (get_lbwmax(ubt) < ubt%lbw+1 .and. ubt%lbw < get_n(ubt)-1)) then
       call set_error(2, info, error); return
    end if
    if (get_lbwmax(wbv) < ubt%lbw .or. get_ubwmax(wbv) < ubt%ubw) then
       call set_error(3, info, error); return
    end if
    if (get_n(ubt) /= get_n(wbv)) then
       call set_error(4, info, error); return
    end if
    call f_d_convert_ubt_to_wbv(ubt%bc, get_n(ubt), ubt%lbw, ubt%ubw, get_lbwmax(ubt), &
         get_ubwmax(ubt), ubt%numrotsu, ubt%jsu, ubt%csu, ubt%ssu, &
         ubt%numrotst, ubt%kst, ubt%cst, ubt%sst, &
         wbv%br, wbv%lbw, wbv%ubw, get_lbwmax(wbv), get_ubwmax(wbv), &
         wbv%numrotsw, wbv%jsw, wbv%csw, wbv%ssw, &
         wbv%numrotsv, wbv%ksv, wbv%csv, wbv%ssv)
    call pop_id(error)
  end subroutine d_convert_ubt_to_wbv

  subroutine f_d_convert_ubt_to_wbv(b_ubt, n, lbw, ubw, lbwmax_ubt, ubwmax_ubt, numrotsu, &
       jsu, csu, ssu, numrotst, kst, cst, sst, &
       b_wbv, lbw_wbv, ubw_wbv, lbwmax_wbv, ubwmax_wbv, &
       numrotsw, jsw, csw, ssw, numrotsv, ksv, csv, ssv)
    real(kind=dp), dimension(lbwmax_ubt+ubwmax_ubt+1,n), intent(inout) :: b_ubt
    integer(kind=int32), intent(in) :: n, lbw, ubw, lbwmax_ubt, ubwmax_ubt, lbwmax_wbv, ubwmax_wbv
    integer(kind=int32), dimension(n), intent(in) :: numrotsu
    integer(kind=int32), dimension(ubwmax_ubt,n), intent(in) :: jsu
    real(kind=dp), dimension(ubwmax_ubt,n), intent(in) :: csu, ssu
    integer(kind=int32), dimension(n), intent(in) :: numrotst
    integer(kind=int32), dimension(n,lbwmax_ubt), intent(in) :: kst
    real(kind=dp), dimension(n,lbwmax_ubt), intent(in) :: cst, sst

    real(kind=dp), dimension(n,lbwmax_wbv+ubwmax_wbv+1), intent(out) :: b_wbv
    integer(kind=int32), dimension(n), intent(out) :: numrotsw
    integer(kind=int32), dimension(lbwmax_wbv,n), intent(out) :: jsw
    real(kind=dp), dimension(lbwmax_wbv,n), intent(out) :: csw, ssw
    integer(kind=int32), dimension(n), intent(out) :: numrotsv
    integer(kind=int32), dimension(n,ubwmax_wbv), intent(out) :: ksv
    real(kind=dp), dimension(n,ubwmax_wbv), intent(out) :: csv, ssv
    integer(kind=int32), intent(out) :: lbw_wbv, ubw_wbv

    integer(kind=int32) :: j, k, k0, k1, ubw1, lbw1, j0, j1
    type(d_rotation) :: rot
    logical :: full_ubw

    b_wbv(:,1:lbw+ubw+1)=0.0_dp; numrotsv=0
    ssv(:,1:ubw)=0.0_dp; csv(:,1:ubw)=0.0_dp
    ksv(:,1:ubw)=0
    numrotsw=0
    ssw(1:lbw,:)=0.0_dp; csw(1:lbw,:)=0.0_dp
    jsw(1:lbw,:)=0
    lbw_wbv=lbw; ubw_wbv=ubw
    if (n == 1) then
       b_wbv(1,1)=b_ubt(1,1);
       lbw_wbv=0; ubw_wbv=0; return
    end if
    ! must allow for temporary fill-in
    if (ubw < n-1) then
       ubw1=ubw+1
       b_ubt=eoshift(b_ubt,-1,dim=1_int64)
       full_ubw=.false.
    else
       ubw1=ubw
       full_ubw=.true.
    end if
    if (lbw < n-1) then
       lbw1=lbw+1
       b_ubt(lbw1+ubw1+1,:)=0.0_dp
    else
       lbw1=lbw
    end if
    ! Do the lower triangular part first
    ! k is the size of the leading principal submatrix
    do k=n-1,2,-1
       ! Apply T_k
       do j=1,numrotst(k)
          call f_tbc_times_rotation(b_ubt,n,lbw1,ubw1,k,0,cst(k,j),-sst(k,j),kst(k,j))
       end do
       ! Rows in which nonzeros have been introduced into the extra subdiagonal.
       j0=max(k+1,lbw+2)
       j1=min(k+lbw,n)
       ! Apply W_{k-1}
       numrotsw(k-1)=max(j1-j0+1,0)
       do j=j0,j1
          rot=lgivens(get_el_bc(b_ubt,ubw1,j-1,j-lbw1),get_el_bc(b_ubt,ubw1,j,j-lbw1))
          jsw(j1-j+1,k-1)=j-1
          csw(j1-j+1,k-1)=rot%cosine; ssw(j1-j+1,k-1)=rot%sine
          call rotation_times_tbc(trp_rot(rot),b_ubt,n,lbw1,ubw1,0,n-k+1,j-1)
       end do
    end do
    ! Upper triangular part
    do k=1,n-2 ! size of trailing principal submatrix
       do j=1,numrotsu(n-k)
          call f_rotation_times_tbc(csu(j,n-k),ssu(j,n-k),b_ubt,n,lbw1,ubw1,&
               n-k,0,jsu(j,n-k))
       end do
       k0=max(n-k+1,ubw+2)
       k1=min(n-k+ubw,n)
       numrotsv(n-(k+1))=max(k1-k0+1,0)
       do j=k0,k1
          rot=rgivens(get_el_bc(b_ubt,ubw1,j-ubw1,j-1), get_el_bc(b_ubt,ubw1,j-ubw1,j))
          ksv(n-(k+1),k1-j+1)=j-1
          csv(n-(k+1),k1-j+1)=rot%cosine; ssv(n-(k+1),k1-j+1)=rot%sine
          call tbc_times_rotation(b_ubt,n,lbw1,ubw1,0,k+1,rot,j-1)
       end do
    end do
    ! Store the results in b_wbv
    if (.not. full_ubw) then
       b_ubt=eoshift(b_ubt,1,dim=1_int64)
    end if
    call bc_to_br(b_ubt, b_wbv, lbw, ubw)
  end subroutine f_d_convert_ubt_to_wbv

  function z_wbv_of_ubt(ubt,error) result(wbv)
    type(z_wbv) :: wbv
    type(z_ubt), intent(in) :: ubt
    type(error_info), intent(inout), optional :: error

    type(z_ubt), allocatable :: ubt1
    type(routine_info), parameter :: info=info_z_wbv_of_ubt

    if (failure(error)) return
    call push_id(info, error)

    wbv=z_new_wbv(get_n(ubt), ubt%lbw, ubt%ubw)
    ubt1=z_new_ubt(get_n(ubt),ubt%lbw+1,  ubt%ubw+1)
    call copy(ubt1,ubt)
    call z_convert_ubt_to_wbv(ubt1,wbv,error)

    call pop_id(error)
  end function z_wbv_of_ubt

  ! Errors:
  ! 0: no error
  ! 1: n<1
  ! 2: Insufficient storage in ubt
  ! 3: Insufficient stroage in wbv
  ! 4: ubt%n /= wbv%n

  subroutine z_convert_ubt_to_wbv(ubt, wbv, error)
    type(z_ubt) :: ubt
    type(z_wbv) :: wbv
    type(error_info), intent(inout), optional :: error
    type(routine_info), parameter :: info=info_z_convert_ubt_to_wbv

    if (failure(error)) return
    call push_id(info, error)

    if (get_n(ubt) < 1) then
       call set_error(1, info, error); return
    end if
    if ((get_ubwmax(ubt) < ubt%ubw+1 .and. ubt%ubw < get_n(ubt)-1) .or. &
         (get_lbwmax(ubt) < ubt%lbw+1 .and. ubt%lbw < get_n(ubt)-1)) then
       call set_error(2, info, error); return
    end if
    if (get_lbwmax(wbv) < ubt%lbw .or. get_ubwmax(wbv) < ubt%ubw) then
       call set_error(3, info, error); return
    end if
    if (get_n(ubt) /= get_n(wbv)) then
       call set_error(4, info, error); return
    end if
    call f_z_convert_ubt_to_wbv(ubt%bc, get_n(ubt), ubt%lbw, ubt%ubw, get_lbwmax(ubt), &
         get_ubwmax(ubt), ubt%numrotsu, ubt%jsu, ubt%csu, ubt%ssu, &
         ubt%numrotst, ubt%kst, ubt%cst, ubt%sst, &
         wbv%br, wbv%lbw, wbv%ubw, get_lbwmax(wbv), get_ubwmax(wbv), &
         wbv%numrotsw, wbv%jsw, wbv%csw, wbv%ssw, &
         wbv%numrotsv, wbv%ksv, wbv%csv, wbv%ssv)
    call pop_id(error)
  end subroutine z_convert_ubt_to_wbv

  subroutine f_z_convert_ubt_to_wbv(b_ubt, n, lbw, ubw, lbwmax_ubt, ubwmax_ubt, numrotsu, &
       jsu, csu, ssu, numrotst, kst, cst, sst, &
       b_wbv, lbw_wbv, ubw_wbv, lbwmax_wbv, ubwmax_wbv, &
       numrotsw, jsw, csw, ssw, numrotsv, ksv, csv, ssv)
    complex(kind=dp), dimension(lbwmax_ubt+ubwmax_ubt+1,n), intent(inout) :: b_ubt
    integer(kind=int32), intent(in) :: n, lbw, ubw, lbwmax_ubt, ubwmax_ubt, lbwmax_wbv, ubwmax_wbv
    integer(kind=int32), dimension(n), intent(in) :: numrotsu
    integer(kind=int32), dimension(ubwmax_ubt,n), intent(in) :: jsu
    real(kind=dp), dimension(ubwmax_ubt,n), intent(in) :: csu
    complex(kind=dp), dimension(ubwmax_ubt,n), intent(in) :: ssu
    integer(kind=int32), dimension(n), intent(in) :: numrotst
    integer(kind=int32), dimension(n,lbwmax_ubt), intent(in) :: kst
    real(kind=dp), dimension(n,lbwmax_ubt), intent(in) :: cst
    complex(kind=dp), dimension(n,lbwmax_ubt), intent(in) :: sst

    complex(kind=dp), dimension(n,lbwmax_wbv+ubwmax_wbv+1), intent(out) :: b_wbv
    integer(kind=int32), dimension(n), intent(out) :: numrotsw
    integer(kind=int32), dimension(lbwmax_wbv,n), intent(out) :: jsw
    real(kind=dp), dimension(lbwmax_wbv,n), intent(out) :: csw
    complex(kind=dp), dimension(lbwmax_wbv,n), intent(out) :: ssw
    integer(kind=int32), dimension(n), intent(out) :: numrotsv
    integer(kind=int32), dimension(n,ubwmax_wbv), intent(out) :: ksv
    real(kind=dp), dimension(n,ubwmax_wbv), intent(out) :: csv
    complex(kind=dp), dimension(n,ubwmax_wbv), intent(out) :: ssv
    integer(kind=int32), intent(out) :: lbw_wbv, ubw_wbv

    integer(kind=int32) :: j, k, k0, k1, ubw1, lbw1, j0, j1
    type(z_rotation) :: rot
    logical :: full_ubw

    b_wbv(:,1:lbw+ubw+1)=(0.0_dp, 0.0_dp); numrotsv=0
    ssv(:,1:ubw)=(0.0_dp, 0.0_dp); csv(:,1:ubw)=0.0_dp
    ksv(:,1:ubw)=0
    numrotsw=0
    ssw(1:lbw,:)=(0.0_dp, 0.0_dp); csw(1:lbw,:)=0.0_dp
    jsw(1:lbw,:)=0
    lbw_wbv=lbw; ubw_wbv=ubw
    if (n == 1) then
       b_wbv(1,1)=b_ubt(1,1);
       lbw_wbv=0; ubw_wbv=0; return
    end if
    ! must allow for temporary fill-in
    if (ubw < n-1) then
       ubw1=ubw+1
       b_ubt=eoshift(b_ubt,-1,dim=1_int64)
       full_ubw=.false.
    else
       ubw1=ubw
       full_ubw=.true.
    end if
    if (lbw < n-1) then
       lbw1=lbw+1
       b_ubt(lbw1+ubw1+1,:)=(0.0_dp, 0.0_dp)
    else
       lbw1=lbw
    end if
    ! Do the lower triangular part first
    ! k is the size of the leading principal submatrix
    do k=n-1,2,-1
       ! Apply T_k
       do j=1,numrotst(k)
          call f_tbc_times_rotation(b_ubt,n,lbw1,ubw1,k,0,cst(k,j),-sst(k,j),kst(k,j))
       end do
       ! Rows in which nonzeros have been introduced into the extra subdiagonal.
       j0=max(k+1,lbw+2)
       j1=min(k+lbw,n)
       ! Apply W_{k-1}
       numrotsw(k-1)=max(j1-j0+1,0)
       do j=j0,j1
          rot=lgivens(get_el_bc(b_ubt,ubw1,j-1,j-lbw1),get_el_bc(b_ubt,ubw1,j,j-lbw1))
          jsw(j1-j+1,k-1)=j-1
          csw(j1-j+1,k-1)=rot%cosine; ssw(j1-j+1,k-1)=rot%sine
          call rotation_times_tbc(trp_rot(rot),b_ubt,n,lbw1,ubw1,0,n-k+1,j-1)
       end do
    end do
    ! Upper triangular part
    do k=1,n-2 ! size of trailing principal submatrix
       do j=1,numrotsu(n-k)
          call f_rotation_times_tbc(csu(j,n-k),ssu(j,n-k),b_ubt,n,lbw1,ubw1,n-k,0,&
               jsu(j,n-k))
       end do
       k0=max(n-k+1,ubw+2)
       k1=min(n-k+ubw,n)
       numrotsv(n-(k+1))=max(k1-k0+1,0)
       do j=k0,k1
          rot=rgivens(get_el_bc(b_ubt,ubw1,j-ubw1,j-1), get_el_bc(b_ubt,ubw1,j-ubw1,j))
          ksv(n-(k+1),k1-j+1)=j-1
          csv(n-(k+1),k1-j+1)=rot%cosine; ssv(n-(k+1),k1-j+1)=rot%sine
          call tbc_times_rotation(b_ubt,n,lbw1,ubw1,0,k+1,rot,j-1)
       end do
    end do
    ! Store the results in b_wbv
    if (.not. full_ubw) then
       b_ubt=eoshift(b_ubt,1,dim=1_int64)
    end if
    call bc_to_br(b_ubt, b_wbv, lbw, ubw)
  end subroutine f_z_convert_ubt_to_wbv

end module mod_convert_ubt_to_wbv

max_stars_repo_name: m-a-stewart/BandGivensWeight

max_issues_repo_name: m-a-stewart/BandGivensWeight

max_forks_repo_name: m-a-stewart/BandGivensWeight


text:
20051030 16:44:05 nbsp Hey Emily! How about helping me make the Outdoor Adventures page good? You know much more about the rafting stuff than I do. :) Users/JackHaskel

20051107 21:41:35 nbsp If you dont vote, youd better not complain. True that. I hate when people complain and they have the right to vote but dont. I cant vote yet, so I still can complain! Users/JohnDudek

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!   EVB-QMDFF - RPMD molecular dynamics and rate constant calculations on
!               black-box generated potential energy surfaces
!
!   Copyright (c) 2021 by Julien Steffen (steffen@pctc.uni-kiel.de)
!                         Stefan Grimme (grimme@thch.uni-bonn.de) (QMDFF code)
!
!   Permission is hereby granted, free of charge, to any person obtaining a
!   copy of this software and associated documentation files (the "Software"),
!   to deal in the Software without restriction, including without limitation
!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
!   and/or sell copies of the Software, and to permit persons to whom the
!   Software is furnished to do so, subject to the following conditions:
!
!   The above copyright notice and this permission notice shall be included in
!   all copies or substantial portions of the Software.
!
!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
!   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
!   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
!   THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
!   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
!   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
!   DEALINGS IN THE SOFTWARE.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!     subroutine blckmgs: Modified Gram-Schmidt orthonormalization
!               of a real matrix
!
!     part of QMDFF
!
subroutine blckmgs(m,n,ndim,darray)
implicit none

integer,intent(in)::m,n,ndim
real(kind=8),dimension(ndim,n),intent(out)::darray
integer::ii,jj,kk,ll,ibsize,nblcks,istrt,jstrt,iend,ncol,ierr
real(kind=8)::tmp
real(kind=8),dimension(:,:),allocatable::smat
real(kind=8)::thr
real(kind=8)::ddot

thr = epsilon(1.0d0) !Threshold for zero vectors
!
!     Block size optimized for Athlon 1200 MHz with 2.0GB memory for
!     matrices up to 5000x5000
ibsize = 60
!
!     Allocate overlap matrix
!
allocate(smat(ibsize,ibsize),stat=ierr)
if(ierr /= 0)  stop 'Mamory allocation error in blckmgs'
!
!     Calculate the number of blocks
!
nblcks = (n+ibsize-1)/ibsize
ibsize = min(n,ibsize)
!
!     Orthogonalize the first block using modified schmidt
!
do ii=1,ibsize

   tmp = ddot(m,darray(1,ii),1,darray(1,ii),1)
!
!     Linear dependence
!
   if (tmp < thr) then
      darray(1:m,ii) = 0.0d0
      cycle
   end if

   tmp = 1.0d0/sqrt(tmp)
   call dscal(m,tmp,darray(1,ii),1)

   do jj=ii+1,ibsize
      tmp = ddot(m,darray(1,ii),1,darray(1,jj),1)
      call daxpy(m,-tmp,darray(1,ii),1,darray(1,jj),1)
   end do
end do
!
!     Loop over remaining blocks
!
do ii=1,nblcks-1
!
!     Initial and final column and number of columns in the block ii+1
!
   istrt = ii*ibsize+1
   iend  = (ii+1)*ibsize
   iend  = min(n,iend)
   ncol  = iend - istrt + 1
!
!     Orthogonalize the block ii+1 against the previous ones
!
   do jj=1,ii
!
!     Initial index of the block jj
!
      jstrt = (jj-1)*ibsize+1
!
!     multiply the matrices/vectors
!
      call dgemm('t','n',ibsize,ncol,m,1.0d0,darray(1,jstrt),ndim,darray(1,istrt),ndim,0.0d0,smat,ibsize)
      call dgemm('n','n',m,ncol,ibsize,-1.0d0,darray(1,jstrt),ndim,smat,ibsize,1.0d0,darray(1,istrt),ndim)
   end do
!
!     Othogonalize vectors on the block ii+1 among themself
!       using modified schmidt
!
   do kk=istrt,iend
      tmp = ddot(m,darray(1,kk),1,darray(1,kk),1)
!
!     Linear dependence
!
      if(tmp < thr) then
         darray(1:m,kk) = 0.0d0
         cycle
      end if
      tmp = 1.0d0/sqrt(tmp)
      call dscal(m,tmp,darray(1,kk),1)
      do ll=kk+1,iend
         tmp = ddot(m,darray(1,kk),1,darray(1,ll),1)
         call daxpy(m,-tmp,darray(1,kk),1,darray(1,ll),1)
      end do
   end do
end do
!
!     Clean up
!
deallocate(smat,stat=ierr)
if(ierr /= 0) stop 'Mamory deallocation error in blckmgs'

return
end subroutine blckmgs

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: chrinide/EVB-QMDFF

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: chrinide/EVB-QMDFF

max_forks_repo_licenses: MIT Unlicense


text:
passive component Telemetry {

  telemetry port tlmOut

  time get port timeGetOut

  @ Channel C1
  telemetry C1: U32
  telemetry C2: U32 format "{}"
  telemetry C3: U32 update on change
  telemetry C4: string
  telemetry C5: I32 \
    low { red -3, orange -2, yellow -1 } \
    high { red 3, orange 2, yellow 1 }

}

ext: fpp

max_stars_repo_name: kevin-f-ortega/fpp

max_issues_repo_name: kevin-f-ortega/fpp

max_forks_repo_name: kevin-f-ortega/fpp


text:
function call(u1,u2) result(s)
character(kind=1,len=*) :: u1,u2
character(kind=1,len=len(u1)+len(u2)) :: s
end function

max_stars_repo_name: xevolver/xevxml

max_issues_repo_name: xevolver/xevxml

max_forks_repo_name: xevolver/xevxml


text:
      SUBROUTINE SLSQPB (M, MEQ, LA, N, X, XL, XU, F, C, G, A, ACC,
     *                   ITER, MODE, R, L, X0, MU, S, U, V, W, IW,
     *                   ALPHA, F0, GS, H1, H2, H3, H4, T, T0, TOL,
     *                   IEXACT, INCONS, IRESET, ITERMX, LINE, N1,
     *                   N2, N3)

C   NONLINEAR PROGRAMMING BY SOLVING SEQUENTIALLY QUADRATIC PROGRAMS

C        -  L1 - LINE SEARCH,  POSITIVE DEFINITE  BFGS UPDATE  -

C                      BODY SUBROUTINE FOR SLSQP

      INTEGER          IW(*), I, IEXACT, INCONS, IRESET, ITER, ITERMX,
     *                 K, J, LA, LINE, M, MEQ, MODE, N, N1, N2, N3

      DOUBLE PRECISION A(LA,N+1), C(LA), G(N+1), L((N+1)*(N+2)/2),
     *                 MU(LA), R(M+N+N+2), S(N+1), U(N+1), V(N+1), W(*),
     *                 X(N), XL(N), XU(N), X0(N),
     *                 DDOT, DNRM2, LINMIN,
     *                 ACC, ALFMIN, ALPHA, F, F0, GS, H1, H2, H3, H4,
     *                 HUN, ONE, T, T0, TEN, TOL, TWO, ZERO

c     dim(W) =         N1*(N1+1) + MEQ*(N1+1) + MINEQ*(N1+1)  for LSQ
c                     +(N1-MEQ+1)*(MINEQ+2) + 2*MINEQ
c                     +(N1+MINEQ)*(N1-MEQ) + 2*MEQ + N1       for LSEI
c                      with MINEQ = M - MEQ + 2*N1  &  N1 = N+1

      DATA             ZERO /0.0D0/, ONE /1.0D0/, ALFMIN /1.0D-1/,
     *                 HUN /1.0D+2/, TEN /1.0D+1/, TWO /2.0D0/

      IF (MODE) 260, 100, 220

  100 ITERMX = ITER
      IF (ACC.GE.ZERO) THEN
          IEXACT = 0
      ELSE
          IEXACT = 1
      ENDIF
      ACC = ABS(ACC)
      TOL = TEN*ACC
      ITER = 0
      IRESET = 0
      N1 = N + 1
      N2 = N1*N/2
      N3 = N2 + 1
      S(1) = ZERO
      MU(1) = ZERO
      CALL DCOPY(N, S(1),  0, S,  1)
      CALL DCOPY(M, MU(1), 0, MU, 1)

C   RESET BFGS MATRIX

  110 IRESET = IRESET + 1
      IF (IRESET.GT.5) GO TO 255
      L(1) = ZERO
      CALL DCOPY(N2, L(1), 0, L, 1)
      J = 1
      DO 120 I=1,N
         L(J) = ONE
         J = J + N1 - I
  120 CONTINUE

C   MAIN ITERATION : SEARCH DIRECTION, STEPLENGTH, LDL'-UPDATE

  130 ITER = ITER + 1
      MODE = 9
      IF (ITER.GT.ITERMX) GO TO 330

C   SEARCH DIRECTION AS SOLUTION OF QP - SUBPROBLEM

      CALL DCOPY(N, XL, 1, U, 1)
      CALL DCOPY(N, XU, 1, V, 1)
      CALL DAXPY(N, -ONE, X, 1, U, 1)
      CALL DAXPY(N, -ONE, X, 1, V, 1)
      H4 = ONE
      CALL LSQ (M, MEQ, N , N3, LA, L, G, A, C, U, V, S, R, W, IW, MODE)

C   AUGMENTED PROBLEM FOR INCONSISTENT LINEARIZATION

      IF (MODE.EQ.6) THEN
          IF (N.EQ.MEQ) THEN
              MODE = 4
          ENDIF
      ENDIF
      IF (MODE.EQ.4) THEN
          DO 140 J=1,M
             IF (J.LE.MEQ) THEN
                 A(J,N1) = -C(J)
             ELSE
                 A(J,N1) = MAX(-C(J),ZERO)
             ENDIF
  140     CONTINUE
          S(1) = ZERO
          CALL DCOPY(N, S(1), 0, S, 1)
          H3 = ZERO
          G(N1) = ZERO
          L(N3) = HUN
          S(N1) = ONE
          U(N1) = ZERO
          V(N1) = ONE
          INCONS = 0
  150     CALL LSQ (M, MEQ, N1, N3, LA, L, G, A, C, U, V, S, R,
     *              W, IW, MODE)
          H4 = ONE - S(N1)
          IF (MODE.EQ.4) THEN
              L(N3) = TEN*L(N3)
              INCONS = INCONS + 1
              IF (INCONS.GT.5) GO TO 330
              GOTO 150
          ELSE IF (MODE.NE.1) THEN
              GOTO 330
          ENDIF
      ELSE IF (MODE.NE.1) THEN
          GOTO 330
      ENDIF

C   UPDATE MULTIPLIERS FOR L1-TEST

      DO 160 I=1,N
         V(I) = G(I) - DDOT(M,A(1,I),1,R,1)
  160 CONTINUE
      F0 = F
      CALL DCOPY(N, X, 1, X0, 1)
      GS = DDOT(N, G, 1, S, 1)
      H1 = ABS(GS)
      H2 = ZERO
      DO 170 J=1,M
         IF (J.LE.MEQ) THEN
             H3 = C(J)
         ELSE
             H3 = ZERO
         ENDIF
         H2 = H2 + MAX(-C(J),H3)
         H3 = ABS(R(J))
         MU(J) = MAX(H3,(MU(J)+H3)/TWO)
         H1 = H1 + H3*ABS(C(J))
  170 CONTINUE

C   CHECK CONVERGENCE

      MODE = 0
      IF (H1.LT.ACC .AND. H2.LT.ACC) GO TO 330
      H1 = ZERO
      DO 180 J=1,M
         IF (J.LE.MEQ) THEN
             H3 = C(J)
         ELSE
             H3 = ZERO
         ENDIF
         H1 = H1 + MU(J)*MAX(-C(J),H3)
  180 CONTINUE
      T0 = F + H1
      H3 = GS - H1*H4
      MODE = 8
      IF (H3.GE.ZERO) GO TO 110

C   LINE SEARCH WITH AN L1-TESTFUNCTION

      LINE = 0
      ALPHA = ONE
      IF (IEXACT.EQ.1) GOTO 210

C   INEXACT LINESEARCH

  190     LINE = LINE + 1
          H3 = ALPHA*H3
          CALL DSCAL(N, ALPHA, S, 1)
          CALL DCOPY(N, X0, 1, X, 1)
          CALL DAXPY(N, ONE, S, 1, X, 1)
          MODE = 1
          GO TO 330
  200         IF (H1.LE.H3/TEN .OR. LINE.GT.10) GO TO 240
              ALPHA = MAX(H3/(TWO*(H3-H1)),ALFMIN)
              GO TO 190

C   EXACT LINESEARCH

  210 IF (LINE.NE.3) THEN
          ALPHA = LINMIN(LINE,ALFMIN,ONE,T,TOL)
          CALL DCOPY(N, X0, 1, X, 1)
          CALL DAXPY(N, ALPHA, S, 1, X, 1)
          MODE = 1
          GOTO 330
      ENDIF
      CALL DSCAL(N, ALPHA, S, 1)
      GOTO 240

C   CALL FUNCTIONS AT CURRENT X

  220     T = F
          DO 230 J=1,M
             IF (J.LE.MEQ) THEN
                 H1 = C(J)
             ELSE
                 H1 = ZERO
             ENDIF
             T = T + MU(J)*MAX(-C(J),H1)
  230     CONTINUE
          H1 = T - T0
          GOTO (200, 210) IEXACT+1

C   CHECK CONVERGENCE

  240 H3 = ZERO
      DO 250 J=1,M
         IF (J.LE.MEQ) THEN
             H1 = C(J)
         ELSE
             H1 = ZERO
         ENDIF
         H3 = H3 + MAX(-C(J),H1)
  250 CONTINUE
      IF ((ABS(F-F0).LT.ACC .OR. DNRM2(N,S,1).LT.ACC) .AND. H3.LT.ACC)
     *   THEN
            MODE = 0
         ELSE
            MODE = -1
         ENDIF
      GO TO 330

C   CHECK relaxed CONVERGENCE in case of positive directional derivative

  255 CONTINUE
      IF ((ABS(F-F0).LT.TOL .OR. DNRM2(N,S,1).LT.TOL) .AND. H3.LT.TOL)
     *   THEN
            MODE = 0
         ELSE
            MODE = 8
         ENDIF
      GO TO 330

C   CALL JACOBIAN AT CURRENT X

C   UPDATE CHOLESKY-FACTORS OF HESSIAN MATRIX BY MODIFIED BFGS FORMULA

  260 DO 270 I=1,N
         U(I) = G(I) - DDOT(M,A(1,I),1,R,1) - V(I)
  270 CONTINUE

C   L'*S

      K = 0
      DO 290 I=1,N
         H1 = ZERO
         K = K + 1
         DO 280 J=I+1,N
            K = K + 1
            H1 = H1 + L(K)*S(J)
  280    CONTINUE
         V(I) = S(I) + H1
  290 CONTINUE

C   D*L'*S

      K = 1
      DO 300 I=1,N
         V(I) = L(K)*V(I)
         K = K + N1 - I
  300 CONTINUE

C   L*D*L'*S

      DO 320 I=N,1,-1
         H1 = ZERO
         K = I
         DO 310 J=1,I - 1
            H1 = H1 + L(K)*V(J)
            K = K + N - J
  310    CONTINUE
         V(I) = V(I) + H1
  320 CONTINUE

      H1 = DDOT(N,S,1,U,1)
      H2 = DDOT(N,S,1,V,1)
      H3 = 0.2D0*H2
      IF (H1.LT.H3) THEN
          H4 = (H2-H3)/(H2-H1)
          H1 = H3
          CALL DSCAL(N, H4, U, 1)
          CALL DAXPY(N, ONE-H4, V, 1, U, 1)
      ENDIF
      CALL LDL(N, L, U, +ONE/H1, V)
      CALL LDL(N, L, V, -ONE/H2, U)

C   END OF MAIN ITERATION

      GO TO 130

C   END OF SLSQPB

  330 END

max_stars_repo_name: robfalck/pyoptsparse

max_issues_repo_name: robfalck/pyoptsparse

max_forks_repo_name: robfalck/pyoptsparse


Program final
	Implicit none
	real*8 :: infty, x, dx, e1, e2, e
	real*8, dimension(2) :: y
	real*8, external :: f_x, g_x
	character (len=10) :: filename
	integer :: n

! parameters
	infty = 4.d0 ! where psi(infty) is considered to be tol

	dx = 0.1d-2/7.d0

	! e1 = -5.3228961107477559
	! e2 = -5.3228961107462656
	! n = 1

	! e1 = -5.3034634104460840
	! e2 = -5.3034634104460832
	! n = 2

	e = -5.3034634104460840
	n = 2

	x = -infty

	y(1) = epsilon(0.d0) ! initial psi at x=infty
	y(2) = epsilon(0.d0) ! initial psi' at x=infty

	write(filename , '("res",i1,".dat")') n
	open(10, file=filename,action="write",status="replace")

	do while( x.le.infty )
		call rk4(x, dx, f_x, g_x, y, e)
		write (10,*) x, y(1)
		x = x + dx
	enddo

	close(10)
end program final

! 4th-order Runge-Kutta subroutine
subroutine rk4(x, dx, df, dg, y, e)
	implicit none
	real*8, external :: df, dg
	real*8, intent(in) :: x, dx, e
	real*8, intent(inout), dimension(2) :: y
	real*8 :: h, k0, k1, k2, k3, l0, l1, l2, l3

	h=dx/2.0D0

	k0 = dx * df(y(1),y(2))
	l0 = dx * dg(y(1),y(2),x,e)
	k1 = dx * df(y(1)+h, y(2)+0.5d0*l0)
	l1 = dx * dg(y(1)+0.5d0*k0,y(2)+0.5d0*l0,x,e)
	k2 = dx * df(y(1)+0.5d0*k1,y(2)+0.5d0*l1)
	l2 = dx * dg(y(1)+0.5d0*k1,y(2)+0.5d0*l1,x,e)
	k3 = dx * df(y(1)+k2,y(2)+l2)
	l3 = dx * dg(y(1)+k2,y(2)+l2,x,e)
	y(1) = y(1) + (k0+2*k1+2*k2+k3)/6.0d0
	y(2) = y(2) + (l0+2*l1+2*l2+l3)/6.0d0
	Return
End subroutine rk4

! function which returns the derivatives (RHS)
real*8 function f_x(y, yp)
! dpsi/dx = psi'
	Implicit none
	real*8 ,intent(in) :: y, yp
	f_x = yp
	Return
End function f_x

real*8 function g_x(y, yp, x, e)
! d^2y/dx^2 = 2(V(x)-e)*y = (4x^4-16x^2-2e)y
	implicit none
	real*8, intent(in) :: y, yp, x, e
!	g_x = ( -4.0d0 - e) * y
	g_x = ( 4.0d0 * x**4.0d0 - 16.0d0 * x**2.0d0 - 2.0d0*e) * y
!	g_x = 2.d0 * ( (x+sqrt(2.d0))**2.d0-8 - e) * y
	Return
End function g_x


text:
      INTEGER FUNCTION IDAMAX(N,DX,INCX)
C
C     FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      DOUBLE PRECISION DX(1),DMAX
      INTEGER I,INCX,IX,N
C
      IDAMAX = 0
      IF( N .LT. 1 ) RETURN
      IDAMAX = 1
      IF(N.EQ.1)RETURN
      IF(INCX.EQ.1)GO TO 20
C
C        CODE FOR INCREMENT NOT EQUAL TO 1
C
      IX = 1
      DMAX = DABS(DX(1))
      IX = IX + INCX
      DO 10 I = 2,N
         IF(DABS(DX(IX)).LE.DMAX) GO TO 5
         IDAMAX = I
         DMAX = DABS(DX(IX))
    5    IX = IX + INCX
   10 CONTINUE
      RETURN
C
C        CODE FOR INCREMENT EQUAL TO 1
C
   20 DMAX = DABS(DX(1))
      DO 30 I = 2,N
         IF(DABS(DX(I)).LE.DMAX) GO TO 30
         IDAMAX = I
         DMAX = DABS(DX(I))
   30 CONTINUE
      RETURN
      END

ext: for


text:
c ***************************************************************
c *                                                             *
c * domain_compute - drive execution of element routines to     *
c *                  compute j and i-integrals for a single     *
c *                  domain                                     *
c *                                                             *
c *                  -> drive MPI to process all elements       *
c *                                                             *
c *                  last update:  5/31/2018 rhd                *
c *                                                             *
c ***************************************************************
c
c
      subroutine dicmj
c
      use global_data, only : myid, iout=>out, nelblk, elblks,
     a                        worker_processor, numprocs, num_threads
c
      use j_data, only : debug_driver, ring_count, symmetric_domain,
     a  print_elem_values, q_element_maps, q_values,
     b  front_element_list, comput_i, e_front, nu_front, comput_j,
     d  cf_traction_flags, expanded_front_nodes, front_q_area,
     e  face_loading, static_j,
     f  domain_min_j, domain_max_j,
     g  domain_min_i, domain_max_i, domain_avg_i, static_min,
     h  static_max, static_avg, ks_from_j, nowring, domain_avg_j,
     i  e33_front, domain_id, j_from_ks, print_totals, j_storage,
     j  i_storage, size_j_values, size_i_values,
     k  j_geonl_formulation, j_linear_formulation,
     l  temperatures_on_model
c
      implicit none
c
c          local declarations.
c
      integer :: now_thread
      integer, external ::  omp_get_thread_num
c
      double precision ::
     &   jtotal, symm_factor, di_value_j, static_di,
     &   static_total, Ki, T11, T33, T13, sign
c
      double precision :: diterms(size_i_values), itotal(size_i_values),
     &   iiterms(size_i_values,size_i_values),
     &   di_value_i(size_j_values)
c
      double precision, allocatable :: diterms_thrds(:,:),
     &                                 iiterms_thrds(:,:,:)
c
      integer :: i, j, skipped_killed, ierr, blk
      integer, save :: message_count=0, message_count2=0
      logical ::  chk_killed
      double precision, parameter :: zero=0.0d0, one=1.0d0,
     &                               two=2.0d0, toler_static_j=1.0d-5
c
c          if running parallel MPI:
c             alert all workers to enter this routine, then send all
c             data required for J-I calculations that is only on root.
c
      call wmpi_alert_slaves( 27 )
      call wmpi_send_jint
c
      if( debug_driver ) write(iout,*) ' >>> entered element driver'
c
c          keep track of ring_count for storing and retrieving output
c
      ring_count = ring_count + 1
c
c          initialize totals for domain and set symmetry factor.
c          print necessary column headers. use arrays to reduce
c          per thread contibutions for J and II terms. saves
c          lots of atomic updates inside block driver
c          routine.
c
      diterms          = zero ! all
      iiterms          = zero !  "
      skipped_killed   = 0
      message_count    = 0
      message_count2   = 0
      symm_factor      = one
      if( symmetric_domain ) symm_factor = two
      if( myid .eq. 0 ) then
         if( print_elem_values ) write(iout,9100)
         if( print_elem_values ) write(iout,9110)
      end if
c
      allocate( diterms_thrds(size_j_values,num_threads),
     &   iiterms_thrds(size_i_values,size_i_values,num_threads) )
      diterms_thrds = zero
      iiterms_thrds = zero
c
c          loop over all element blocks. inside block checks
c          if element is in current domain
c
c$OMP PARALLEL DO PRIVATE( now_thread, blk ) ! all else shared
      do blk = 1, nelblk
         if( elblks(2,blk) .ne. myid ) cycle
         now_thread = omp_get_thread_num() + 1
         call dicmj_do_a_blk( blk, now_thread,
     &      diterms_thrds(1,now_thread), iiterms_thrds(1,1,now_thread),
     &      message_count, message_count2  )
      end do
c$OMP END PARALLEL DO
c
c          reduce all thread contributions to J and II terms
c
      do i = 1, num_threads
        diterms = diterms + diterms_thrds(:,i)
        iiterms = iiterms + iiterms_thrds(:,:,i)
      end do
c
      call wmpi_reduce_vec( diterms, size_j_values )
      call wmpi_reduce_vec( iiterms(1,1), size_i_values**2 )
      call wmpi_redint( skipped_killed )
      call wmpi_redlog( face_loading )
c
      if( worker_processor ) then
        if( allocated( q_values ) ) deallocate( q_values )
        if( allocated( q_element_maps ) ) deallocate( q_element_maps )
        if( allocated( front_element_list ) )
     &      deallocate( front_element_list )
        if( allocated( expanded_front_nodes ) )
     &      deallocate( expanded_front_nodes )
        return
      end if
c
c             done with this domain computation. print values
c             as required based on user specified flags.
c
c             only rank 0 runs this code. see deallocate/return
c             above for workers
c
c             J-integral results:
c
      if( comput_j ) call dicmj_finalize_j_terms
c
c             interaction-integral results:
c             calculate K from relationship with interaction integral
c
c             the 8 terms for I-integral results are stored
c             in array e_iresults(8,7) as follows:
c
c                               value    auxiliary field
c
c                  iiterms(i,1): KI       plane stress
c                  iiterms(i,2): KI       plane stress
c                  iiterms(i,3): KII      plane stress
c                  iiterms(i,4): KII      plane stress
c                  iiterms(i,5): KIII     anti-plane shear
c                  iiterms(i,6): T11,T33  plane stress
c                  iiterms(i,7): T11,T33  plane strain
c                  iiterms(i,8): T13      anti-plane strain
c
      if( comput_i ) call dicmj_finalize_i_terms
c
c             compute J-values from K-values...
c
c             including KI plane stress, KII plane stress, and KIII
c
      j_from_ks(ring_count,1) = i_storage(ring_count,11,1)**2 / e_front
     &                        + i_storage(ring_count,11,3)**2 / e_front
     &                        + i_storage(ring_count,11,5)**2
     &                        * (one + nu_front ) / e_front
c
c             including KI plane strain, KII plane strain, and KIII
c
      j_from_ks(ring_count,2) = i_storage(ring_count,11,2)**2
     &                          * (one - nu_front**2 ) / e_front
     &                        + i_storage(ring_count,11,4)**2
     &                          * (one - nu_front**2 ) / e_front
     &                        + i_storage(ring_count,11,5)**2
     &                        * (one + nu_front ) / e_front
c
c             print sum of values from all elements in domain
c
      call dicmj_print_values
c
      return
c
 9100 format(/,1x,25x,'domain integral components',
     &  /,      1x,25x,'--------------------------')
 9110 format(7x,'element',5x,'dm1',9x,'dm2',9x,'dm3',9x,'dm4',
     &    9x,'dm5',9x,'dm6',9x,'dm7',9x,'dm8' )
c
      contains
c     ========
c
c***************************************************************
c                                                              *
c      dicmj -> contains dicmj_finalize_j_terms                *
c                                                              *
c        written by: rhd                                       *
c        last modified: 3/30/2018 rhd                          *
c                                                              *
c***************************************************************
c
      subroutine dicmj_finalize_j_terms
      implicit none
c
      jtotal = diterms(1) + diterms(2) + diterms(3) + diterms(4) +
     &         diterms(5) + diterms(6) + diterms(7) + diterms(8)
      if( abs(diterms(3)+diterms(4)).lt.toler_static_j*abs(jtotal) )
     &     static_j = .true.
      static_total = diterms(1) + diterms(2) + diterms(5)
     &             + diterms(6) + diterms(7) + diterms(8)
      di_value_j   = jtotal / front_q_area
      static_di    = static_total / front_q_area
      domain_min_j = min( domain_min_j, di_value_j )
      domain_max_j = max( domain_max_j, di_value_j )
      domain_avg_j = domain_avg_j + di_value_j
      static_min   = min( static_min, static_di  )
      static_max   = max( static_max, static_di  )
      static_avg   = static_avg + static_di
c
c             calculate stress intensity factor K, from J, for:
c
c               1. plane stress where loading is pure mode I or mode II.
c               2. plane strain where loading is pure mode I or mode II.
c               3. anti-plane shear where loading is pure mode III.
c
c             if J is negative, calculate K from |J|, and then change sign.
c
      sign = di_value_j / abs(di_value_j)
c
      ks_from_j(ring_count,1) = sqrt(abs(di_value_j) * e_front)*sign
c
      ks_from_j(ring_count,2) =
     &     sqrt(abs(di_value_j) * e_front/(one - nu_front**2))*sign
c
      ks_from_j(ring_count,3) =
     &     sqrt(abs(di_value_j) * e_front/(one + nu_front)) * sign
c
c          store J-integral values for output in didriv.f
c
      j_storage(ring_count,1)   = dble( nowring )
      j_storage(ring_count,2:9) = diterms(1:8)
      j_storage(ring_count,10)  = di_value_j
      j_storage(ring_count,11)  = dble(skipped_killed)
c
      return
      end subroutine dicmj_finalize_j_terms
c
c***************************************************************
c                                                              *
c      dicmj -> contains dicmj_finalize_i_terms                *
c                                                              *
c        written by: rhd                                       *
c        last modified: 3/30/2018 rhd                          *
c                                                              *
c***************************************************************
c
      subroutine dicmj_finalize_i_terms
      implicit none
c
      do j = 1, 8
         itotal(j) = iiterms(1,j) + iiterms(2,j) + iiterms(3,j)
     &             + iiterms(4,j) + iiterms(5,j) + iiterms(6,j)
     &             + iiterms(7,j) + iiterms(8,j)
         di_value_i(j) = itotal(j) / front_q_area
c
c          calculate stress intensity factors and T-stress
c
         Ki  = zero
         T11 = zero
         T33 = zero
         T13 = zero
c
c          KI (plane stress auxiliary fields)
c
         if( j .eq. 1 ) Ki = di_value_i(j) * e_front / two
c
c          KI (plane strain auxiliary fields)
c
         if( j .eq. 2 ) Ki = di_value_i(j) * e_front
     &                     / ( two * (one - nu_front**2) )
c
c          KII (plane stress auxiliary fields)
c
         if( j .eq. 3 ) Ki = di_value_i(j) * e_front / two
c
c          KII (plane strain auxiliary fields)
c
         if( j .eq. 4 ) Ki = di_value_i(j) * e_front
     &                     / ( two * (one - nu_front**2) )
c
c          KIII (anti-plane shear auxiliary fields)
c
         if( j .eq. 5 ) Ki = di_value_i(j) * e_front
     &                     / (two * (one + nu_front))
c
c          T11, T33 (plane stress auxiliary fields)
c
         if( j .eq. 6 ) then
            T11 = e_front * di_value_i(j)
            Ki  = T11
            T33 = zero
            domain_min_i(9) = T33
            domain_max_i(9) = T33
            domain_avg_i(9) = T33
         end if
c
c          T11, T33 (plane strain auxiliary fields)
c
         if( j .eq. 7 ) then
           T11 = e_front / (one - nu_front**2 )
     &         * ( di_value_i(j) + nu_front * e33_front )
           T33 = e33_front * e_front + nu_front * T11
           Ki  = T11
           domain_min_i(10) = min( domain_min_i(10), T33 )
           domain_max_i(10) = max( domain_max_i(10), T33 )
           domain_avg_i(10) = domain_avg_i(10) + T33
         end if
c
c             T13 (anti-plane shear auxiliary fields)
c
         if( j .eq. 8 ) then
            T13 = di_value_i(j) * e_front / (two * (one + nu_front ))
            Ki  = T13
         end if
c
c             determine maximum and minimum values
c
         domain_min_i(j) = min( domain_min_i(j), Ki )
         domain_max_i(j) = max( domain_max_i(j), Ki )
         domain_avg_i(j) = domain_avg_i(j) + Ki
c
c             store i-integral values for output in didriv.f
c
         i_storage(ring_count,1,j)   = dble( nowring )
         i_storage(ring_count,2:9,j) = iiterms(1:8,j)
         i_storage(ring_count,10,j)  = di_value_i(j)
         i_storage(ring_count,11,j)  = Ki
         if( j.eq.7 ) i_storage(ring_count,12,j) = T33
         i_storage(ring_count,13,j)  = dble(skipped_killed)
      end do  !  on i
c
      return
      end subroutine dicmj_finalize_i_terms
c***************************************************************
c                                                              *
c      dicmj -> contains dicmj_print_values                    *
c                                                              *
c        written by: rhd                                       *
c        last modified: 3/30/2018 rhd                          *
c                                                              *
c***************************************************************
c
      subroutine dicmj_print_values
      implicit none

      if ( print_elem_values ) then
         write(iout,9160)
         if( comput_j ) write(iout,9170) "J     ", (diterms(i),i=1,8),
     &        jtotal
         if( comput_i ) then
            write(iout,9170) "I_KI  ", (iiterms(i,1),i=1,8), itotal(1)
            write(iout,9170) "I_KI  ", (iiterms(i,2),i=1,8), itotal(2)
            write(iout,9170) "I_KII ", (iiterms(i,3),i=1,8), itotal(3)
            write(iout,9170) "I_KII ", (iiterms(i,4),i=1,8), itotal(4)
            write(iout,9170) "I_KIII", (iiterms(i,5),i=1,8), itotal(5)
            write(iout,9170) "I_T11 ", (iiterms(i,6),i=1,8), itotal(6)
            write(iout,9170) "I_T11 ", (iiterms(i,7),i=1,8), itotal(7)
            write(iout,9170) "I_T13 ", (iiterms(i,8),i=1,8), itotal(8)
         end if
      end if
c
c             print sum of values from all elements in domain,
c             and J, I for domain
c
      if( print_totals ) then   ! should not be used now (3/9/2018, rhd)
         write(iout,9180)
         if( nowring .eq. 0 ) then
            if( comput_j .and. comput_i ) write(iout,9185)
            if( comput_j .and. .not. comput_i ) write(iout,9186)
            if( comput_i .and. .not. comput_j ) write(iout,9187)
            if( comput_j ) then
               write(iout,9220) "J     ", domain_id(1:8),
     &              (diterms(j),j=1,8), di_value_j, skipped_killed
            end if
            if( comput_i ) then
               write(iout,9220) "I_KI  ", domain_id(1:8),
     &              (iiterms(i,1),i=1,8), di_value_i(1), skipped_killed
               write(iout,9220) "I_KI  ", domain_id(1:8),
     &              (iiterms(i,2),i=1,8), di_value_i(2), skipped_killed
               write(iout,9220) "I_KII ", domain_id(1:8),
     &              (iiterms(i,3),i=1,8), di_value_i(3), skipped_killed
               write(iout,9220) "I_KII ", domain_id(1:8),
     &              (iiterms(i,4),i=1,8), di_value_i(4), skipped_killed
               write(iout,9220) "I_KIII", domain_id(1:8),
     &              (iiterms(i,5),i=1,8), di_value_i(5), skipped_killed
               write(iout,9220) "I_T11 ", domain_id(1:8),
     &              (iiterms(i,6),i=1,8), di_value_i(6), skipped_killed
               write(iout,9220) "I_T11 ", domain_id(1:8),
     &              (iiterms(i,7),i=1,8), di_value_i(7), skipped_killed
               write(iout,9220) "I_T13 ", domain_id(1:8),
     &              (iiterms(i,8),i=1,8), di_value_i(8), skipped_killed
            end if
         end if
c
         if( nowring.gt.0 ) then
            if( comput_j .and.       comput_i ) write(iout,9250)
            if( comput_j .and. .not. comput_i ) write(iout,9270)
            if( comput_i .and. .not. comput_j ) write(iout,9290)
            if( comput_j ) then
               write(iout,9300) "J     ", nowring, (diterms(i),i=1,8),
     &              di_value_j, skipped_killed
            end if
            if( comput_i ) then
               write(iout,9300) "I_KI  ", nowring,
     &              (iiterms(i,1),i=1,8), di_value_i(1), skipped_killed
               write(iout,9300) "I_KI  ", nowring,
     &              (iiterms(i,2),i=1,8), di_value_i(2), skipped_killed
               write(iout,9300) "I_KII ", nowring,
     &              (iiterms(i,3),i=1,8), di_value_i(3), skipped_killed
               write(iout,9300) "I_KII ", nowring,
     &              (iiterms(i,4),i=1,8), di_value_i(4), skipped_killed
               write(iout,9300) "I_KIII", nowring,
     &              (iiterms(i,5),i=1,8), di_value_i(5), skipped_killed
               write(iout,9300) "I_T11 ", nowring,
     &              (iiterms(i,6),i=1,8), di_value_i(6), skipped_killed
               write(iout,9300) "I_T11 ", nowring,
     &              (iiterms(i,7),i=1,8), di_value_i(7), skipped_killed
               write(iout,9300) "I_T13 ", nowring,
     &              (iiterms(i,8),i=1,8), di_value_i(8), skipped_killed
            end if
         end if
      end if
c
      return
c
 9100 format(/,1x,25x,'domain integral components',
     &  /,      1x,25x,'--------------------------')
 9110 format(7x,'element',5x,'dm1',9x,'dm2',9x,'dm3',9x,'dm4',
     &    9x,'dm5',9x,'dm6',9x,'dm7',9x,'dm8' )
 9130 format(/,5x,' gauss point strains for element ',i7,', gpt ',i2,
     &       ':',/, 3(10x,3(e11.4,2x),/))
 9140 format(1x,a,1x,i7,8(1x,e11.4))
 9150 format('')
 9160 format(/,1x,'element totals:',3x,'dm1',9x,'dm2',9x,'dm3',9x,'dm4',
     &    9x,'dm5',9x,'dm6',9x,'dm7',9x,'dm8',8x,'total')
 9170 format(1x,a,7x,9(1x,e11.4))
 9180 format(/,1x,25x,'domain integral components',
     &  /,      1x,25x,'--------------------------')
 9185 format(8x,'domain',8x,'dm1',9x,'dm2',9x,'dm3',9x,'dm4',
     &    9x,'dm5',9x,'dm6',9x,'dm7',9x,'dm8',5x,'total J,I',
     &    2x,'killed ele' )
 9186 format(8x,'domain',8x,'dm1',9x,'dm2',9x,'dm3',9x,'dm4',
     &    9x,'dm5',9x,'dm6',9x,'dm7',9x,'dm8',7x,'total J',
     &    2x,'killed ele' )
 9187 format(8x,'domain',8x,'dm1',9x,'dm2',9x,'dm3',9x,'dm4',
     &    9x,'dm5',9x,'dm6',9x,'dm7',9x,'dm8',7x,'total I',
     &    2x,'killed ele' )
 9220 format(1x,a,1x,a8,9(1x,e11.4),2x,'(',i3,')')
 9250 format(8x,'domain',6x,'dm1',9x,'dm2',9x,'dm3',9x,'dm4',
     &    9x,'dm5',9x,'dm6',9x,'dm7',9x,'dm8',5x,'total J,I',
     &    2x,'killed ele' )
 9270 format(8x,'domain',6x,'dm1',9x,'dm2',9x,'dm3',9x,'dm4',
     &    9x,'dm5',9x,'dm6',9x,'dm7',9x,'dm8',7x,'total J',
     &    2x,'killed ele' )
 9290 format(8x,'domain',6x,'dm1',9x,'dm2',9x,'dm3',9x,'dm4',
     &    9x,'dm5',9x,'dm6',9x,'dm7',9x,'dm8',7x,'total I',
     &    2x,'killed ele' )
 9300 format(1x,a,1x,i7,9(1x,e11.4),2x,'(',i3,')')
 9400 format(/,1x,'area under q-function along crack front:  ',e11.4,
     &       /,1x,'length along crack front for this domain: ',e11.4,/)
c
      end subroutine dicmj_print_values
c
      end subroutine dicmj

c     ****************************************************************
c     *                                                              *
c     *                      dicmj_do_a_blk                          *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 6/23/2018 rhd              *
c     *                                                              *
c     *      Processing 1 domain and output values. Run MPI          *
c     *      parallel over (mpi) domains, threads over blocks        *
c     *      in each rank. In large models, many (MPI) domains       *
c     *      and element blocks have no elements within the J-II     *
c     *      domain and are skipped                                  *
c     *                                                              *
c     ****************************************************************
c
      subroutine dicmj_do_a_blk( blk, now_thread, diterms,
     &                           iiterms, message_count,
     &                           message_count2 )
c
      use global_data, only : myid, iout=>out, noelem, nonode, bits,
     a                        iprops, props, lprops, elblks, mxndel,
     b                        nstrs, mxgp, mxelpr, nstr,
     c                        worker_processor, numprocs,
     d                        scoord=>c, sdispl=>u, sveloc=>v,
     e                        saccel=>a
      use main_data, only   : incmap, elems_to_blocks, incid,
     a                        fgm_node_values, eq_node_force_indexes,
     b                        trn, trnmat,
     c                        temper_elems, temper_nodes,
     d                        temper_nodes_ref, eq_node_forces,
     e                        eq_node_force_len
      use elem_block_data, only : rot_n1_blocks,
     a                            urcs_n_blocks, cdest_blocks,
     b                            edest_blocks, eps_n_blocks,
     c                            initial_state_data
      use j_data, only : debug_driver, symmetric_domain,
     a  print_elem_values, q_element_maps, q_values, fgm_e, fgm_nu,
     b  block_seg_curves, seg_snode_nu, swd_at_nodes, strain_at_nodes,
     c  ignore_face_loads, omit_crack_front_elems, front_element_list,
     d  comput_i, cf_tractions, domain_rot, debug_elements,
     e  one_point_rule, e_front, nu_front, comput_j,
     f  cf_traction_flags, face_loading, process_temperatures,
     g  snode_alpha_ij, seg_snode_e, displ_grad_at_nodes,
     h  size_j_values, size_i_values, front_list_length,
     i  j_geonl_formulation, j_linear_formulation,
     j  process_initial_state, temperatures_on_model
c
      implicit none
c
c          parameters
c
      integer :: blk, now_thread, message_count, message_count2
      double precision :: diterms(size_j_values),
     &                    iiterms(size_i_values,size_i_values)
c
c          local declarations. includes pointers to simplify
c          addressing into blocks.
c
      double precision ::
     &   e_coord(3,mxndel), e_displ(3,mxndel), e_vel(3,mxndel),
     &   e_accel(3,mxndel), q_element(mxndel),
     &   e_stress(nstrs*mxgp), e_rots(9*mxgp), e_force(3,mxndel),
     &   symm_factor, e_node_temps(mxndel),
     &   elem_uniform_temp, e_alpha_ij(6,mxndel), eq_load_modifier,
     &   dummy, elem_nod_strains(6,mxndel), elem_nod_swd(mxndel),
     &   elem_nod_displ_derivs(9,mxndel), e_strain(9,mxgp), sign,
     &   e(mxndel), nu(mxndel), cf_load(3), e_W_is(mxgp)
c
      double precision :: e_jresults(size_j_values),
     &   e_iresults(size_i_values,size_i_values)
c
      integer :: eps_offset, gpn, node, node_id, etype,
     &           skipped_killed, ierr, elemno, num_enodes,
     &           num_gpts, rel_elem, i, j, i_elem,
     &           offset, element, numrow_sig, span, felem
      integer :: snodes(mxndel)
      integer, dimension (:,:), pointer :: edest, cdest
      real :: e_props(mxelpr)
      logical, external :: dibmck
      logical ::  geonl, chk_killed, omit_J7_J8_front_elems, ok,
     &            elem_temps, face_loads_possible, is_front_element,
     &            front_elem_flag, seg_curves_flag, brick, linearform,
     &            process_swd_derivs, process_strain_derivs,
     &            process_grad_derivs, process_blk
      double precision, parameter :: zero=0.0d0, one=1.0d0,
     &                               two=2.0d0, toler_static_j=1.0d-5
c
      span  = elblks(0,blk)
      felem = elblks(1,blk)
c
c          loop over all elements in this block. is there work to
c          do: (1) element is in domain, and (2) not killed
c
      process_blk = .false.
c
      do i_elem = 1, span
         elemno = felem + i_elem - 1
         if( .not. dibmck( elemno, q_element_maps, bits ) ) cycle
         if( chk_killed( elemno ) ) cycle
         process_blk = .true.
         exit
      end do
      if( .not. process_blk ) return
c
      symm_factor      = one
      if( symmetric_domain ) symm_factor = two
      process_swd_derivs    = size( swd_at_nodes ) == nonode
      process_strain_derivs = size( strain_at_nodes, 2 ) == nonode
      process_grad_derivs   = size( displ_grad_at_nodes, 2 ) == nonode
c
      do i_elem = 1, span
c
c          handle only elements in J-II domain
c
         elemno = felem + i_elem - 1
         if( .not. dibmck( elemno, q_element_maps, bits ) ) cycle
         if( chk_killed( elemno ) ) then
c$OMP ATOMIC UPDATE
            skipped_killed = skipped_killed + 1
            cycle
         end if
c
         etype      = iprops(1,elemno)
         brick      = etype .ge. 1  .and. etype .le. 5
         if( .not. brick ) then
c$OMP ATOMIC UPDATE
             message_count = message_count + 1
             if( message_count <=10 ) then
               write(iout,9200) elemno
               if( message_count == 10 ) write(iout,9210)
             end if
             cycle
         end if
c
         num_enodes = iprops(2,elemno)
         num_gpts   = iprops(6,elemno)
         geonl      = lprops(18,elemno)
         linearform = .not. geonl
         ok = .true.
         if( j_linear_formulation ) then  ! element must match domain
            if( geonl ) ok = .false.
         end if
         if( j_geonl_formulation ) then
            if( linearform ) ok = .false.
         end if
         if( .not. ok ) then
c$OMP ATOMIC UPDATE
             message_count2 = message_count2 + 1
             if( message_count2 <=10 ) then
               write(iout,9230) elemno
               if( message_count2 == 10 ) write(iout,9240)
             end if
             cycle
         end if
c
         rel_elem   = elemno - felem + 1
         cdest      => cdest_blocks(blk)%ptr
         edest      => edest_blocks(blk)%ptr
c
c          build copy of element coordinates, displacements,
c          velocities, accelerations, q-values from global data.
c          zero load terms that correspond to constrained dof.
c          if we have temperature loadings, get temps for element
c          nodes (node values + element uniform value). rotate
c          displacements, velocities and accelerations from
c          constraint-compatible to global coordinates. obtain
c          element nodal values of stress work density and displ
c          gradients structure node-value arrays.
c
         call dicmj_element_setup
c
c          crack-face tractions. the contribution to J arising
c          from crack-face tractions is calculated element-by-element
c          using the equivalent nodal loads on each element face
c          computed during setup of load step for applied pressures,
c          surface tractions and body forces. here we just get
c          the element equiv. nodal forces. they are passed to
c          element J routine which sorts out the loaded face.
c
c          a similar procedure is used to compute the contribution
c          to I.
c
         face_loads_possible = .not. ignore_face_loads .and.
     &                               eq_node_force_len .gt. 0
         if( face_loads_possible ) then
            if( eq_node_force_indexes(elemno) .ne. 0 ) then
               call vec_ops( e_force,  !    copy routine
     &              eq_node_forces(eq_node_force_indexes(elemno)),
     &              dummy, 3*num_enodes, 5 )
            end if
         end if
c
         is_front_element = .false.
         do i = 1, front_list_length
           if( front_element_list(i) == elemno ) then
             is_front_element = .true.
             exit
            end if
         end do
         omit_J7_J8_front_elems = omit_crack_front_elems .and.
     &                            is_front_element
c
c          gather element stresses
c
         offset = (rel_elem-1) * nstrs * num_gpts
         associate( urcs_n =>  urcs_n_blocks(blk)%ptr )
         do i = 1, nstrs * num_gpts
            e_stress(i) = urcs_n(offset+i)
         end do
         end associate

c          gather element properties
c
         e_props(1:mxelpr) = props(1:mxelpr,elemno)
c
c          gather 3x3 rotation matrices from polar decompositions
c          at each gauss point of element for geonl
c
         if( geonl ) then
           offset =  (rel_elem-1) * 9 * num_gpts
           associate( rots => rot_n1_blocks(blk)%ptr )
           do i = 1, 9 * num_gpts
             e_rots(i) = rots(offset+i)
           end do
           end associate
         end if
c
         e_strain = zero ! all terms
c
c          for I computations,
c          gather element strains at integration points.
c          and store in tensor form. convert engineering
c          (total) strains of off-diagonal terms to tensor strains.
c
c          for crack-face loading, make a copy of the tractions
c          input by the user in the domain definition. if none
c          was input, dielem_c.f automatically uses the equivalent
c          nodal loads used to solve the boundary value problem.
c
         if( comput_i ) call dicmj_more_i_setup
c
c          set flag for crack-front element.
c
         front_elem_flag = is_front_element
c
c          call the element routine to compute contribution
c          to the j-integral and or i-integral for domain
c
         element    = elemno  !  protect variables
         numrow_sig = nstrs   !       "
         call dielem( e_coord, q_element, e_displ, e_vel, e_accel,
     a                e_force, e_stress, e_props, e_props, e_props,
     b                e_rots, domain_rot, e_jresults, e_node_temps,
     c                elem_temps, e_alpha_ij, ierr, element, geonl,
     d                numrow_sig, snodes, e_W_is, elem_nod_swd,
     e                elem_nod_strains, elem_nod_displ_derivs,
     f                omit_J7_J8_front_elems, fgm_e,
     f                fgm_nu, e_strain, e, e_front, nu, nu_front,
     g                e_iresults, cf_traction_flags, cf_load,
     h                front_elem_flag, seg_curves_flag )
c
         if( comput_j ) then
            do i = 1, size_j_values
               e_jresults(i) = e_jresults(i) * symm_factor
               diterms(i)    = diterms(i) + e_jresults(i)
            end do
            if( print_elem_values .and. comput_j )
     &        write(iout,9140) "J     ", elemno,
     &                         e_jresults(1:size_j_values)
         end if
c
         if( comput_i ) then
            do i = 1, size_i_values
               do j = 1, size_i_values
                  e_iresults(i,j) = e_iresults(i,j) * symm_factor
                  iiterms(i,j)    = iiterms(i,j) + e_iresults(i,j)
               end do
            end do
         end if
c
         if( comput_i .and. symmetric_domain ) then
            iiterms(1:size_i_values,3:5) = zero
            iiterms(1:size_i_values,8)   = zero
         end if
c
         if( print_elem_values .and. comput_i ) then
            write(iout,9140) "I_KI  ", elemno, (e_iresults(i,1),i=1,8)
            write(iout,9140) "I_KI  ", elemno, (e_iresults(i,2),i=1,8)
            write(iout,9140) "I_KII ", elemno, (e_iresults(i,3),i=1,8)
            write(iout,9140) "I_KII ", elemno, (e_iresults(i,4),i=1,8)
            write(iout,9140) "I_KIII", elemno, (e_iresults(i,5),i=1,8)
            write(iout,9140) "I_T11 ", elemno, (e_iresults(i,6),i=1,8)
            write(iout,9140) "I_T11 ", elemno, (e_iresults(i,7),i=1,8)
            write(iout,9140) "I_T13 ", elemno, (e_iresults(i,8),i=1,8)
            write(iout,9150)
         end if
      end do  ! **** on all elements in block ****

      return
 9140 format(1x,a,1x,i7,8(1x,e11.4))
 9150 format('')
 9200 format(/,'>> Warning: only brick elements supported for domain',
     &  ' integrals. element: ',i8,
     & /,    '           element skipped.')

 9210 format(/,'>> Warning: *** no more messages about bricks ***')
 9230 format(/,'>> Warning: all elements in domain must have the same',
     &  ' formulation (small strain, larger strain. element: ',i8,
     & /,    '           element skipped.')
 9240 format(/,'>> Warning: *** no more messages about formulation',
     &   ' mismatch ***')

      contains
c     ========
c
c***************************************************************
c                                                              *
c      dicmj_do_a_blk -> contains dicmj_element_setup          *
c                                                              *
c        written by: rhd                                       *
c        last modified: 6/23/2018 rhd                          *
c                                                              *
c***************************************************************
c
      subroutine dicmj_element_setup
      implicit none
c
      integer :: incpos, k1, k2, i1, i2, i3, enode, snode
      logical, parameter :: local_debug = .true.
c
      seg_curves_flag = .false.
c
      incpos     = incmap(elemno)
      k1         = num_enodes
      k2         = 2*num_enodes
c
!DIR$ VECTOR ALIGNED
      do enode = 1, num_enodes
         e_coord(1,enode) = scoord(cdest(enode,rel_elem))
         e_coord(2,enode) = scoord(cdest(enode+k1,rel_elem))
         e_coord(3,enode) = scoord(cdest(enode+k2,rel_elem))
         snode            = incid(incpos+enode-1)
         snodes(enode)    = snode
         q_element(enode) = q_values(snode)
         i1               = edest(enode,   rel_elem)
         i2               = edest(k1+enode,rel_elem)
         i3               = edest(k2+enode,rel_elem)
         e_displ(1,enode) = sdispl(i1)
         e_displ(2,enode) = sdispl(i2)
         e_displ(3,enode) = sdispl(i3)
         e_vel(1,enode)   = sveloc(i1)
         e_vel(2,enode)   = sveloc(i2)
         e_vel(3,enode)   = sveloc(i3)
         e_accel(1,enode) = saccel(i1)
         e_accel(2,enode) = saccel(i2)
         e_accel(3,enode) = saccel(i3)
         e(enode)         = props(7,elemno)
         nu(enode)        = props(8,elemno)
         if( fgm_e )  e(enode)  = fgm_node_values(snode,1)
         if( fgm_nu ) nu(enode) = fgm_node_values(snode,2)
         if( temperatures_on_model .and.
     &       allocated( block_seg_curves ) ) then
            if( block_seg_curves(blk) ) then
               seg_curves_flag = .true.
               e(enode)        = seg_snode_e(snode)
               nu(enode)       = seg_snode_nu(snode)
            end if
         end if
c
c         if necessary, transform nodal values from
c         constraint-compatible coordinates to global
c         coordinates.
c
         if( trn(snode) ) then
            call di_trans_nodvals( e_displ(1,enode),
     &                             trnmat(snode)%mat(1,1))
            call di_trans_nodvals( e_vel(1,enode),
     &                             trnmat(snode)%mat(1,1))
            call di_trans_nodvals( e_accel(1,enode),
     &                             trnmat(snode)%mat(1,1))
         end if
c
         e_force(1,enode)          = zero
         e_force(2,enode)          = zero
         e_force(3,enode)          = zero
         e_node_temps(enode)       = zero
         e_alpha_ij(1:6,enode)     = zero
         elem_nod_swd(enode)       = zero
         elem_nod_strains(1:6,enode) = zero
         elem_nod_displ_derivs(1:9,enode) = zero
c
      end do ! on enode
c

      if( process_swd_derivs ) then
!DIR$ VECTOR ALIGNED
        do enode = 1, num_enodes
          snode = incid(incpos+enode-1)
          elem_nod_swd(enode) = swd_at_nodes(snode)
        end do
      end if
c
      if( process_strain_derivs ) then
        do enode = 1, num_enodes
          snode = incid(incpos+enode-1)
!DIR$ VECTOR ALIGNED
          elem_nod_strains(1:6,enode) = strain_at_nodes(1:6,snode)
        end do
      end if
c
      if( process_grad_derivs ) then
        do enode = 1, num_enodes
          snode = incid(incpos+enode-1)
!DIR$ VECTOR ALIGNED
          elem_nod_displ_derivs(1:9,enode) =
     &               displ_grad_at_nodes(1:9,snode)
        end do
      end if
c
c        for problems with temperature loads, pull out
c        current temperatures of element nodes and set the thermal
c        expansion coefficients at element nodes. remove reference
c        temperatures of model nodes. the global temperatures flag
c        refers to increment temps for step just completed. cannot use
c        it since early application of temps may be followed by only
c        mechanical loads
c
      elem_temps = .false.
      if( temperatures_on_model ) then
        elem_uniform_temp = temper_elems(elemno)
!DIR$ VECTOR ALIGNED
        do enode = 1, num_enodes
         snode = incid(incpos + enode - 1)
         e_node_temps(enode) = temper_nodes(snode) + elem_uniform_temp
     &                        - temper_nodes_ref(snode)
         if( abs(e_node_temps(enode)) > zero ) elem_temps = .true.
        end do
        if( allocated( snode_alpha_ij ) ) then
         do enode = 1, num_enodes
          snode = incid( incpos + enode - 1 )
!DIR$ VECTOR ALIGNED
          e_alpha_ij(1:6,enode) = snode_alpha_ij(snode,1:6)
         end do
        end if
      end if
c
c        setup vector of energy density values at element
c        integration points at user-defined initial state
c        if required. otherwise set zero vector to eliminate
c        logic in element computational routines.
c
      if( process_initial_state ) then
        do gpn = 1, num_gpts
         e_W_is(gpn) =
     &     initial_state_data(blk)%W_plastic_nis_block(i_elem,gpn)
        end do

        if( local_debug .and.
     &      (elemno == 8955 .or. elemno == 9060 ) ) then
          write(iout,*) ' @ 1. dicmj_element_setup. elem: ',elemno
          do gpn = 1, num_gpts
           write(iout,9000)  gpn, e_W_is(gpn)
          end do
        end if
      else
          e_W_is = zero
      end if
 9000 format(10x,i3,f15.8)
c
      return
      end subroutine dicmj_element_setup
c
c***************************************************************
c                                                              *
c      dicmj -> contains dicmj_more_i_setup                    *
c                                                              *
c        written by: rhd                                       *
c        last modified: 3/30/2018 rhd                          *
c                                                              *
c***************************************************************
c
      subroutine dicmj_more_i_setup
      implicit none

c
c          gather element strains at integration points.
c          and store in tensor form. convert engineering
c          (total) strains of off-diagonal terms to tensor strains.
c
      eps_offset = (rel_elem - 1) * nstr * num_gpts
      associate( eps => eps_n_blocks(blk)%ptr )
c
!DIR$ VECTOR ALIGNED
      do gpn = 1, num_gpts
         e_strain(1,gpn) = eps(eps_offset + 1)
         e_strain(2,gpn) = eps(eps_offset + 4)/two
         e_strain(3,gpn) = eps(eps_offset + 6)/two
         e_strain(4,gpn) = eps(eps_offset + 4)/two
         e_strain(5,gpn) = eps(eps_offset + 2)
         e_strain(6,gpn) = eps(eps_offset + 5)/two
         e_strain(7,gpn) = eps(eps_offset + 6)/two
         e_strain(8,gpn) = eps(eps_offset + 5)/two
         e_strain(9,gpn) = eps(eps_offset + 3)
         eps_offset      = eps_offset + nstr
      end do
c
      end associate
c
      if( debug_driver ) then
         do gpn = 1, num_gpts
            write(iout,9130) elemno, gpn, (e_strain(i,gpn),i=1,9)
         end do
      end if
c
c          for crack-face loading, make a copy of the tractions
c          input by the user in the domain definition. if none
c          was input, dielem_c.f automatically uses the equivalent
c          nodal loads used to solve the boundary value problem.
c
      cf_load(1) = cf_tractions(1)
      cf_load(2) = cf_tractions(2)
      cf_load(3) = cf_tractions(3)
c
      return
 9130 format(/,5x,' gauss point strains for element ',i7,', gpt ',i2,
     &       ':',/, 3(10x,3(e11.4,2x),/))

      end subroutine  dicmj_more_i_setup
c
      end subroutine dicmj_do_a_blk

c***************************************************************
c                                                              *
c subroutine to transform a 3x1 vector of constraint-          *
c compatible-coordinate-system values to global-coordinate-    *
c system values. the routine premultiplies the incoming        *
c constraint-compatible values by the transpose of the stored  *
c global-to-constraint coordinate system rotation matrix.      *
c                                                              *
c                written by: mcw                               *
c                last modified: 3/26/2018 rhd                  *
c                                                              *
c***************************************************************
c
      subroutine di_trans_nodvals( vec, transmat )
      implicit none
      double precision :: vec(3), transmat(3,3)
      double precision :: tempvec(3)
c
      tempvec(1) = transmat(1,1) * vec(1) + transmat(2,1) * vec(2)
     &           + transmat(3,1) * vec(3)
      tempvec(2) = transmat(1,2) * vec(1) + transmat(2,2) * vec(2)
     &           + transmat(3,2) * vec(3)
      tempvec(3) = transmat(1,3) * vec(1) + transmat(2,3) * vec(2)
     &           + transmat(3,3) * vec(3)
c
      vec(1) = tempvec(1)
      vec(2) = tempvec(2)
      vec(3) = tempvec(3)
c
      return
      end