!> @brief This module contains the MODFLOW 6 BMI
!!
!! This BMI interface matches the CSDMS standard, with a few modifications:
!!
!! - This interface will build into a shared library that can be called from other
!!   executables and scripts, not necessarily written in Fortran. Therefore we have
!!   omitted the type-boundness of the routines, since they cannot have the
!!   bind(C,"...") attribute.
!! - MODFLOW has internal data arrays with rank > 1 that we would like to expose. An
!!   example would be access to data in the BOUND array of GWF boundary packages (BndType).
!!   The get_value_ptr calls below support this, returning a C-style pointer to the arrays
!!   and methods have been added to query the variable's rank and shape.
!!
!! Note on style: BMI apparently uses underscores, we use underscores in some
!! places but camelcase in other. Since this is a dedicated BMI interface module,
!! we'll use underscores here as well.
!<
module mf6bmi
  use mf6bmiUtil
  use mf6bmiError
  use Mf6CoreModule
  use TdisModule, only: kper, kstp
  use iso_c_binding, only: c_int, c_char, c_double, C_NULL_CHAR, c_loc, c_ptr, &
                           c_f_pointer
  use KindModule, only: DP, I4B, LGP
  use ConstantsModule, only: LENMEMPATH, LENVARNAME
  use MemoryManagerModule, only: mem_setptr, get_mem_elem_size, get_isize, &
                                 get_mem_rank, get_mem_shape, get_mem_type, &
                                 memorylist, get_from_memorylist
  use MemoryTypeModule, only: MemoryType
  use MemoryHelperModule, only: create_mem_address
  use SimVariablesModule, only: simstdout, istdout
  use InputOutputModule, only: getunit
  implicit none

  integer(c_int), bind(C, name="ISTDOUTTOFILE") :: istdout_to_file = 1  !< output control: =0 to screen, >0 to file
  !DIR$ ATTRIBUTES DLLEXPORT :: istdout_to_file

contains

  function bmi_get_component_name(name) result(bmi_status) bind(C, name="get_component_name")
    !DIR$ ATTRIBUTES DLLEXPORT :: bmi_get_component_name
    ! -- dummy variables
    character(kind=c_char), intent(out) :: name(BMI_LENCOMPONENTNAME)
    integer(kind=c_int) :: bmi_status !< BMI status code
    ! -- local variables

    name = string_to_char_array('MODFLOW 6', 9)
    bmi_status = BMI_SUCCESS

  end function bmi_get_component_name

  !> @brief Initialize the computational core
  !!
  !! It is required to have the MODFLOW 6 configuration file 'mfsim.nam'
  !! available in the current working directory when calling this routine.
  !!
  !! NOTE: initialization should always be matched with a call to
  !! bmi_finalize(). However, we currently do not support the reinitialization
  !! of a model in the same memory space... You would have to create a new
  !! process for that.
  !<
  function bmi_initialize() result(bmi_status) bind(C, name="initialize")
    !DIR$ ATTRIBUTES DLLEXPORT :: bmi_initialize
    ! -- dummy variables
    integer(kind=c_int) :: bmi_status !< BMI status code
    ! -- local variables

    if (istdout_to_file > 0) then
      ! -- open stdout file mfsim.stdout
      istdout = getunit()
      !
      ! -- set STDOUT to a physical file unit
      open (unit=istdout, file=simstdout)
    end if
    !
    ! -- initialize MODFLOW 6
    call Mf6Initialize()

    bmi_status = BMI_SUCCESS

  end function bmi_initialize

  !> @brief Perform a computational time step
  !!
  !! It will prepare the timestep, call the calculation routine
  !! on all the solution groups in the simulation, and finalize
  !! the timestep by printing out diagnostics and writing output.
  !! It can be called in succession to perform multiple steps up
  !! to the simulation's end time is reached.
  !<
  function bmi_update() result(bmi_status) bind(C, name="update")
    !DIR$ ATTRIBUTES DLLEXPORT :: bmi_update
    ! -- dummy variables
    integer(kind=c_int) :: bmi_status !< BMI status code
    ! -- local variables
    logical :: hasConverged

    hasConverged = Mf6Update()

    bmi_status = BMI_SUCCESS
  end function bmi_update

  !> @brief Clean up the initialized simulation
  !!
  !! Performs teardown tasks for the initialized simulation, this
  !! call should match the call to bmi_initialize()
  !<
  function bmi_finalize() result(bmi_status) bind(C, name="finalize")
    !DIR$ ATTRIBUTES DLLEXPORT :: bmi_finalize
    ! -- modules
    use SimVariablesModule, only: iforcestop
    ! -- dummy variables
    integer(kind=c_int) :: bmi_status !< BMI status code

    ! we don't want a full stop() here, this disables it:
    iforcestop = 0
    call Mf6Finalize()

    bmi_status = BMI_SUCCESS

  end function bmi_finalize

  !> @brief Get the start time of the simulation
  !!
  !! As MODFLOW currently does not have internal time, this will be
  !! returning 0.0 for now. New version...
  !<
  function get_start_time(start_time) result(bmi_status) bind(C, name="get_start_time")
    !DIR$ ATTRIBUTES DLLEXPORT :: get_start_time
    ! -- dummy variables
    real(kind=c_double), intent(out) :: start_time !< start time
    integer(kind=c_int) :: bmi_status           !< BMI status code

    start_time = 0.0_DP
    bmi_status = BMI_SUCCESS

  end function get_start_time

  !> @brief Get the end time of the simulation
  !!
  !! As MODFLOW does currently does not have internal time, this will be
  !! equal to the total runtime.
  !<
  function get_end_time(end_time) result(bmi_status) bind(C, name="get_end_time")
    !DIR$ ATTRIBUTES DLLEXPORT :: get_end_time
    ! -- modules
    use TdisModule, only: totalsimtime
    ! -- dummy variables
    real(kind=c_double), intent(out) :: end_time !< end time
    integer(kind=c_int) :: bmi_status            !< BMI status code

    end_time = totalsimtime
    bmi_status = BMI_SUCCESS

  end function get_end_time

  !> @brief Get the current time of the simulation
  !!
  !! As MODFLOW currently does not have internal time, this will be
  !! equal to the time passed w.r.t. the start time of the simulation.
  !<
  function get_current_time(current_time) result(bmi_status) bind(C, name="get_current_time")
    !DIR$ ATTRIBUTES DLLEXPORT :: get_current_time
    ! -- modules
    use TdisModule, only: totim
    ! -- dummy variables
    real(kind=c_double), intent(out) :: current_time !< current time
    integer(kind=c_int) :: bmi_status             !< BMI status code

    current_time = totim
    bmi_status = BMI_SUCCESS

  end function get_current_time

  !> @brief Get the time step for the simulation
  !!
  !! Note that the returned value may vary between and within stress periods,
  !! depending on your time discretization settings in the TDIS package.
  !<
  function get_time_step(time_step) result(bmi_status) bind(C, name="get_time_step")
    !DIR$ ATTRIBUTES DLLEXPORT :: get_time_step
    ! -- modules
    use TdisModule, only: delt
    ! -- dummy variables
    real(kind=c_double), intent(out) :: time_step  !< current time step
    integer(kind=c_int) :: bmi_status           !< BMI status code

    time_step = delt
    bmi_status = BMI_SUCCESS

  end function get_time_step

  !> @brief Get the number of input variables in the simulation
  !!
  !! This concerns all variables stored in the memory manager
  !<
  function get_input_item_count(count) result(bmi_status) bind(C, name="get_input_item_count")
    !DIR$ ATTRIBUTES DLLEXPORT :: get_input_item_count
    ! -- dummy variables
    integer(kind=c_int), intent(out) :: count !< the number of input variables
    integer(kind=c_int) :: bmi_status         !< BMI status code

    count = memorylist%count()

    bmi_status = BMI_SUCCESS

  end function get_input_item_count

  !> @brief Get the number of output variables in the simulation
  !!
  !! This concerns all variables stored in the memory manager
  !<
  function get_output_item_count(count) result(bmi_status) bind(C, name="get_output_item_count")
    !DIR$ ATTRIBUTES DLLEXPORT :: get_output_item_count
    ! -- dummy variables
    integer(kind=c_int), intent(out) :: count !< the number of output variables
    integer(kind=c_int) :: bmi_status         !< BMI status code

    count = memorylist%count()

    bmi_status = BMI_SUCCESS

  end function get_output_item_count

  !> @brief Returns all input variables in the simulation
  !!
  !! This functions returns the full address for all variables in the
  !! memory manager
  !!
  !! The array @p c_names should be pre-allocated of proper size:
  !!
  !! size = BMI_LENVARADDRESS *  get_input_item_count()
  !!
  !! The strings will be written contiguously with stride equal to
  !! BMI_LENVARADDRESS and nul-terminated where the trailing spaces start:
  !!
  !! c_names = 'variable_address_1\x00 ... variable_address_2\x00 ... ' etc.
  !<
  function get_input_var_names(c_names) result(bmi_status) bind(C, name="get_input_var_names")
    !DIR$ ATTRIBUTES DLLEXPORT :: get_input_var_names
    ! -- dummy variables
    character(kind=c_char, len=1), intent(inout) :: c_names(*) !< array with memory paths for input variables
    integer(kind=c_int) :: bmi_status                         !< BMI status code
    ! -- local variables
    integer(I4B) :: imem, start, i
    type(MemoryType), pointer :: mt => null()
    character(len=LENMEMADDRESS) :: var_address

    start = 1
    do imem = 1, memorylist%count()
      mt => memorylist%Get(imem)
      var_address = create_mem_address(mt%path, mt%name)
      do i = 1, len(trim(var_address))
        c_names(start + i - 1) = var_address(i:i)
      end do
      c_names(start + i) = c_null_char
      start = start + BMI_LENVARADDRESS
    end do

    bmi_status = BMI_SUCCESS

  end function get_input_var_names

  !> @brief Returns all output variables in the simulation
  !!
  !! This function works analogously to get_input_var_names(),
  !! and currently returns the same set of memory variables,
  !! which is all of them!
  !<
  function get_output_var_names(c_names) result(bmi_status) bind(C, name="get_output_var_names")
    !DIR$ ATTRIBUTES DLLEXPORT :: get_output_var_names
    ! -- dummy variables
    character(kind=c_char, len=1), intent(inout) :: c_names(*) !< array with memory paths for output variables
    integer(kind=c_int) :: bmi_status                         !< BMI status code
    ! -- local variables
    integer(I4B) :: imem, start, i
    type(MemoryType), pointer :: mt => null()
    character(len=LENMEMADDRESS) :: var_address

    start = 1
    do imem = 1, memorylist%count()
      mt => memorylist%Get(imem)
      var_address = create_mem_address(mt%path, mt%name)
      do i = 1, len(trim(var_address))
        c_names(start + i - 1) = var_address(i:i)
      end do
      c_names(start + i) = c_null_char
      start = start + BMI_LENVARADDRESS
    end do

    bmi_status = BMI_SUCCESS

  end function get_output_var_names

  !> @brief Get the size (in bytes) of a single element of a variable
  !<
  function get_var_itemsize(c_var_address, var_size) result(bmi_status) bind(C, name="get_var_itemsize")
    !DIR$ ATTRIBUTES DLLEXPORT :: get_var_itemsize
    ! -- dummy variables
    character(kind=c_char), intent(in) :: c_var_address(*) !< memory address string of the variable
    integer(kind=c_int), intent(out) :: var_size                        !< size of the element in bytes
    integer(kind=c_int) :: bmi_status                       !< BMI status code
    ! -- local variables
    character(len=LENMEMPATH) :: mem_path
    character(len=LENVARNAME) :: var_name
    logical(LGP) :: valid

    bmi_status = BMI_SUCCESS

    call split_address(c_var_address, mem_path, var_name, valid)
    if (.not. valid) then
      bmi_status = BMI_FAILURE
      return
    end if

    call get_mem_elem_size(var_name, mem_path, var_size)
    if (var_size == -1) bmi_status = BMI_FAILURE

  end function get_var_itemsize

  !> @brief Get size of the variable, in bytes
  !<
  function get_var_nbytes(c_var_address, var_nbytes) result(bmi_status) bind(C, name="get_var_nbytes")
    !DIR$ ATTRIBUTES DLLEXPORT :: get_var_nbytes
    ! -- dummy variables
    character(kind=c_char), intent(in) :: c_var_address(*) !< memory address string of the variable
    integer(kind=c_int), intent(out) :: var_nbytes                      !< size in bytes
    integer(kind=c_int) :: bmi_status                       !< BMI status code
    ! -- local variables
    integer(I4B) :: var_size, isize
    character(len=LENMEMPATH) :: mem_path
    character(len=LENVARNAME) :: var_name
    logical(LGP) :: valid

    bmi_status = BMI_SUCCESS

    call split_address(c_var_address, mem_path, var_name, valid)
    if (.not. valid) then
      bmi_status = BMI_FAILURE
      return
    end if

    call get_mem_elem_size(var_name, mem_path, var_size)
    if (var_size == -1) bmi_status = BMI_FAILURE
    call get_isize(var_name, mem_path, isize)
    if (isize == -1) bmi_status = BMI_FAILURE

    var_nbytes = var_size*isize

  end function get_var_nbytes

  !> @brief Copy the double precision values of a variable into the array
  !!
  !! The copied variable us located at @p c_var_address. The caller should
  !! provide @p c_arr_ptr pointing to an array of the proper shape (the
  !! BMI function get_var_shape() can be used to create it). Multi-dimensional
  !! arrays are supported.
  !<
  function get_value_double(c_var_address, c_arr_ptr) result(bmi_status) bind(C, name="get_value_double")
    !DIR$ ATTRIBUTES DLLEXPORT :: get_value_double
    ! -- modules
    use MemorySetHandlerModule, only: on_memory_set
    ! -- dummy variables
    character(kind=c_char), intent(in) :: c_var_address(*) !< memory address string of the variable
    type(c_ptr), intent(in) :: c_arr_ptr                    !< pointer to the double precision array
    integer(kind=c_int) :: bmi_status                                   !< BMI status code
    ! -- local variables
    character(len=LENMEMPATH) :: mem_path
    character(len=LENVARNAME) :: var_name
    logical(LGP) :: valid
    integer(I4B) :: rank
    real(DP), pointer :: src_ptr, tgt_ptr
    real(DP), dimension(:), pointer, contiguous :: src1D_ptr, tgt1D_ptr
    real(DP), dimension(:, :), pointer, contiguous :: src2D_ptr, tgt2D_ptr
    real(DP), dimension(:, :, :), pointer, contiguous :: src3D_ptr, tgt3D_ptr
    integer(I4B) :: i, j, k

    bmi_status = BMI_SUCCESS

    call split_address(c_var_address, mem_path, var_name, valid)
    if (.not. valid) then
      bmi_status = BMI_FAILURE
      return
    end if

    ! convert pointer and copy data from memory manager into
    ! the passed array, using loops to avoid stack overflow
    rank = -1
    call get_mem_rank(var_name, mem_path, rank)

    if (rank == 0) then
      call mem_setptr(src_ptr, var_name, mem_path)
      call c_f_pointer(c_arr_ptr, tgt_ptr)
      tgt_ptr = src_ptr
    else if (rank == 1) then
      call mem_setptr(src1D_ptr, var_name, mem_path)
      call c_f_pointer(c_arr_ptr, tgt1D_ptr, shape(src1D_ptr))
      do i = 1, size(tgt1D_ptr)
        tgt1D_ptr(i) = src1D_ptr(i)
      end do
    else if (rank == 2) then
      call mem_setptr(src2D_ptr, var_name, mem_path)
      call c_f_pointer(c_arr_ptr, tgt2D_ptr, shape(src2D_ptr))
      do j = 1, size(tgt2D_ptr, 2)
        do i = 1, size(tgt2D_ptr, 1)
          tgt2D_ptr(i, j) = src2D_ptr(i, j)
        end do
      end do
    else if (rank == 3) then
      call mem_setptr(src3D_ptr, var_name, mem_path)
      call c_f_pointer(c_arr_ptr, tgt3D_ptr, shape(src3D_ptr))
      do k = 1, size(tgt3D_ptr, 3)
        do j = 1, size(tgt3D_ptr, 2)
          do i = 1, size(tgt3D_ptr, 1)
            tgt3D_ptr(i, j, k) = src3D_ptr(i, j, k)
          end do
        end do
      end do
    else
      write (bmi_last_error, fmt_unsupported_rank) trim(var_name)
      call report_bmi_error(bmi_last_error)
      bmi_status = BMI_FAILURE
      return
    end if

  end function get_value_double

  !> @brief Copy the integer values of a variable into the array
  !!
  !! The copied variable us located at @p c_var_address. The caller should
  !! provide @p c_arr_ptr pointing to an array of the proper shape (the
  !! BMI function get_var_shape() can be used to create it). Multi-dimensional
  !! arrays are supported.
  !<
  function get_value_int(c_var_address, c_arr_ptr) result(bmi_status) bind(C, name="get_value_int")
    !DIR$ ATTRIBUTES DLLEXPORT :: get_value_int
    ! -- modules
    use MemorySetHandlerModule, only: on_memory_set
    ! -- dummy variables
    character(kind=c_char), intent(in) :: c_var_address(*) !< memory address string of the variable
    type(c_ptr), intent(in) :: c_arr_ptr                    !< pointer to the integer array
    integer(kind=c_int) :: bmi_status                                   !< BMI status code
    ! -- local variables
    character(len=LENMEMPATH) :: mem_path
    character(len=LENVARNAME) :: var_name
    logical(LGP) :: valid
    integer(I4B) :: rank
    integer(I4B), pointer :: src_ptr, tgt_ptr
    integer(I4B), dimension(:), pointer, contiguous :: src1D_ptr, tgt1D_ptr
    integer(I4B), dimension(:, :), pointer, contiguous :: src2D_ptr, tgt2D_ptr
    integer(I4B), dimension(:, :, :), pointer, contiguous :: src3D_ptr, tgt3D_ptr
    integer(I4B) :: i, j, k

    bmi_status = BMI_SUCCESS

    call split_address(c_var_address, mem_path, var_name, valid)
    if (.not. valid) then
      bmi_status = BMI_FAILURE
      return
    end if

    ! convert pointer and copy data from memory manager into
    ! the passed array, using loops to avoid stack overflow
    rank = -1
    call get_mem_rank(var_name, mem_path, rank)

    if (rank == 0) then
      call mem_setptr(src_ptr, var_name, mem_path)
      call c_f_pointer(c_arr_ptr, tgt_ptr)
      tgt_ptr = src_ptr
    else if (rank == 1) then
      call mem_setptr(src1D_ptr, var_name, mem_path)
      call c_f_pointer(c_arr_ptr, tgt1D_ptr, shape(src1D_ptr))
      do i = 1, size(tgt1D_ptr)
        tgt1D_ptr(i) = src1D_ptr(i)
      end do
    else if (rank == 2) then
      call mem_setptr(src2D_ptr, var_name, mem_path)
      call c_f_pointer(c_arr_ptr, tgt2D_ptr, shape(src2D_ptr))
      do j = 1, size(tgt2D_ptr, 2)
        do i = 1, size(tgt2D_ptr, 1)
          tgt2D_ptr(i, j) = src2D_ptr(i, j)
        end do
      end do
    else if (rank == 3) then
      call mem_setptr(src3D_ptr, var_name, mem_path)
      call c_f_pointer(c_arr_ptr, tgt3D_ptr, shape(src3D_ptr))
      do k = 1, size(tgt3D_ptr, 3)
        do j = 1, size(tgt3D_ptr, 2)
          do i = 1, size(tgt3D_ptr, 1)
            tgt3D_ptr(i, j, k) = src3D_ptr(i, j, k)
          end do
        end do
      end do
    else
      write (bmi_last_error, fmt_unsupported_rank) trim(var_name)
      call report_bmi_error(bmi_last_error)
      bmi_status = BMI_FAILURE
      return
    end if

  end function get_value_int

  !> @brief Get a pointer to the array of double precision numbers
  !!
  !! The array is located at @p c_var_address. There is no copying of data involved.
  !! Multi-dimensional arrays are supported and the get_var_rank() function
  !! can be used to get the variable's dimensionality, and get_var_shape() for
  !! its shape.
  !<
  function get_value_ptr_double(c_var_address, c_arr_ptr) result(bmi_status) bind(C, name="get_value_ptr_double")
    !DIR$ ATTRIBUTES DLLEXPORT :: get_value_ptr_double
    ! -- dummy variables
    character(kind=c_char), intent(in) :: c_var_address(*) !< memory address string of the variable
    type(c_ptr), intent(inout) :: c_arr_ptr                 !< pointer to the array
    integer(kind=c_int) :: bmi_status                       !< BMI status code
    ! -- local variables
    character(len=LENMEMPATH) :: mem_path
    character(len=LENVARNAME) :: var_name
    logical(LGP) :: valid
    real(DP), pointer :: scalar_ptr
    real(DP), dimension(:), pointer, contiguous :: array_ptr
    real(DP), dimension(:, :), pointer, contiguous :: array2D_ptr
    real(DP), dimension(:, :, :), pointer, contiguous :: array3D_ptr
    integer(I4B) :: rank

    bmi_status = BMI_SUCCESS

    call split_address(c_var_address, mem_path, var_name, valid)
    if (.not. valid) then
      bmi_status = BMI_FAILURE
      return
    end if

    rank = -1
    call get_mem_rank(var_name, mem_path, rank)
    if (rank == 0) then
      call mem_setptr(scalar_ptr, var_name, mem_path)
      c_arr_ptr = c_loc(scalar_ptr)
    else if (rank == 1) then
      call mem_setptr(array_ptr, var_name, mem_path)
      c_arr_ptr = c_loc(array_ptr)
    else if (rank == 2) then
      call mem_setptr(array2D_ptr, var_name, mem_path)
      c_arr_ptr = c_loc(array2D_ptr)
    else if (rank == 3) then
      call mem_setptr(array3D_ptr, var_name, mem_path)
      c_arr_ptr = c_loc(array3D_ptr)
    else
      write (bmi_last_error, fmt_unsupported_rank) trim(var_name)
      call report_bmi_error(bmi_last_error)
      bmi_status = BMI_FAILURE
      return
    end if

  end function get_value_ptr_double

  !> @brief Get a pointer to the array of integer numbers
  !!
  !! The array is located at @p c_var_address. There is no copying of data involved.
  !! Multi-dimensional arrays are supported and the get_var_rank() function
  !! can be used to get the variable's dimensionality.
  !<
  function get_value_ptr_int(c_var_address, c_arr_ptr) result(bmi_status) bind(C, name="get_value_ptr_int")
    !DIR$ ATTRIBUTES DLLEXPORT :: get_value_ptr_int
    ! -- dummy variables
    character(kind=c_char), intent(in) :: c_var_address(*) !< memory address string of the variable
    type(c_ptr), intent(inout) :: c_arr_ptr                 !< pointer to the array
    integer(kind=c_int) :: bmi_status                       !< BMI status code
    ! -- local variables
    character(len=LENMEMPATH) :: mem_path
    character(len=LENVARNAME) :: var_name
    logical(LGP) :: valid
    integer(I4B) :: rank
    integer(I4B), pointer :: scalar_ptr
    integer(I4B), dimension(:), pointer, contiguous :: array_ptr
    integer(I4B), dimension(:, :), pointer, contiguous :: array2D_ptr
    integer(I4B), dimension(:, :, :), pointer, contiguous :: array3D_ptr

    bmi_status = BMI_SUCCESS

    call split_address(c_var_address, mem_path, var_name, valid)
    if (.not. valid) then
      bmi_status = BMI_FAILURE
      return
    end if

    rank = -1
    call get_mem_rank(var_name, mem_path, rank)

    if (rank == 0) then
      call mem_setptr(scalar_ptr, var_name, mem_path)
      c_arr_ptr = c_loc(scalar_ptr)
    else if (rank == 1) then
      call mem_setptr(array_ptr, var_name, mem_path)
      c_arr_ptr = c_loc(array_ptr)
    else if (rank == 2) then
      call mem_setptr(array2D_ptr, var_name, mem_path)
      c_arr_ptr = c_loc(array2D_ptr)
    else if (rank == 3) then
      call mem_setptr(array3D_ptr, var_name, mem_path)
      c_arr_ptr = c_loc(array3D_ptr)
    else
      write (bmi_last_error, fmt_unsupported_rank) trim(var_name)
      call report_bmi_error(bmi_last_error)
      bmi_status = BMI_FAILURE
      return
    end if

  end function get_value_ptr_int

  !> @brief Set new values for a variable of type double
  !!
  !! The array pointed to by @p c_arr_ptr can have rank equal to 0, 1, or 2
  !! and should have a C-style layout, which is particularly important for
  !! rank > 1.
  !<
  function set_value_double(c_var_address, c_arr_ptr) result(bmi_status) bind(C, name="set_value_double")
    !DIR$ ATTRIBUTES DLLEXPORT :: set_value_double
    ! -- modules
    use MemorySetHandlerModule, only: on_memory_set
    ! -- dummy variables
    character(kind=c_char), intent(in) :: c_var_address(*) !< memory address string of the variable
    type(c_ptr), intent(in) :: c_arr_ptr                   !< pointer to the double precision array
    integer(kind=c_int) :: bmi_status                      !< BMI status code
    ! -- local variables
    character(len=LENMEMPATH) :: mem_path
    character(len=LENVARNAME) :: var_name
    logical(LGP) :: valid
    integer(I4B) :: rank
    real(DP), pointer :: src_ptr, tgt_ptr
    real(DP), dimension(:), pointer, contiguous :: src1D_ptr, tgt1D_ptr
    real(DP), dimension(:, :), pointer, contiguous :: src2D_ptr, tgt2D_ptr
    integer(I4B) :: i, j
    integer(I4B) :: status

    bmi_status = BMI_SUCCESS

    call split_address(c_var_address, mem_path, var_name, valid)
    if (.not. valid) then
      bmi_status = BMI_FAILURE
      return
    end if

    ! convert pointer and copy, using loops to avoid stack overflow
    rank = -1
    call get_mem_rank(var_name, mem_path, rank)

    if (rank == 0) then
      call mem_setptr(tgt_ptr, var_name, mem_path)
      call c_f_pointer(c_arr_ptr, src_ptr)
      tgt_ptr = src_ptr
    else if (rank == 1) then
      call mem_setptr(tgt1D_ptr, var_name, mem_path)
      call c_f_pointer(c_arr_ptr, src1D_ptr, shape(tgt1D_ptr))
      do i = 1, size(tgt1D_ptr)
        tgt1D_ptr(i) = src1D_ptr(i)
      end do
    else if (rank == 2) then
      call mem_setptr(tgt2D_ptr, var_name, mem_path)
      call c_f_pointer(c_arr_ptr, src2D_ptr, shape(tgt2D_ptr))
      do j = 1, size(tgt2D_ptr, 2)
        do i = 1, size(tgt2D_ptr, 1)
          tgt2D_ptr(i, j) = src2D_ptr(i, j)
        end do
      end do
    else
      write (bmi_last_error, fmt_unsupported_rank) trim(var_name)
      call report_bmi_error(bmi_last_error)
      bmi_status = BMI_FAILURE
      return
    end if

    ! trigger event:
    call on_memory_set(var_name, mem_path, status)
    if (status /= 0) then
      ! something went terribly wrong here, aborting
      write (bmi_last_error, fmt_invalid_mem_access) trim(var_name)
      call report_bmi_error(bmi_last_error)
      bmi_status = BMI_FAILURE
      return
    end if

  end function set_value_double

  !> @brief Set new values for a variable of type integer
  !!
  !! The array pointed to by @p c_arr_ptr can have rank equal to 0, 1, or 2.
  !<
  function set_value_int(c_var_address, c_arr_ptr) result(bmi_status) bind(C, name="set_value_int")
    !DIR$ ATTRIBUTES DLLEXPORT :: set_value_int
    ! -- modules
    use MemorySetHandlerModule, only: on_memory_set
    ! -- dummy variables
    character(kind=c_char), intent(in) :: c_var_address(*) !< memory address string of the variable
    type(c_ptr), intent(in) :: c_arr_ptr                   !< pointer to the integer array
    integer(kind=c_int) :: bmi_status                      !< BMI status code
    ! -- local variables
    character(len=LENMEMPATH) :: mem_path
    character(len=LENVARNAME) :: var_name
    logical(LGP) :: valid
    integer(I4B) :: rank
    integer(I4B), pointer :: src_ptr, tgt_ptr
    integer(I4B), dimension(:), pointer, contiguous :: src1D_ptr, tgt1D_ptr
    integer(I4B), dimension(:, :), pointer, contiguous :: src2D_ptr, tgt2D_ptr
    integer(I4B) :: i, j
    integer(I4B) :: status

    bmi_status = BMI_SUCCESS

    call split_address(c_var_address, mem_path, var_name, valid)
    if (.not. valid) then
      bmi_status = BMI_FAILURE
      return
    end if

    ! convert pointer and copy, using loops to avoid stack overflow
    rank = -1
    call get_mem_rank(var_name, mem_path, rank)

    if (rank == 0) then
      call mem_setptr(tgt_ptr, var_name, mem_path)
      call c_f_pointer(c_arr_ptr, src_ptr)
      tgt_ptr = src_ptr
    else if (rank == 1) then
      call mem_setptr(tgt1D_ptr, var_name, mem_path)
      call c_f_pointer(c_arr_ptr, src1D_ptr, shape(tgt1D_ptr))
      do i = 1, size(tgt1D_ptr)
        tgt1D_ptr(i) = src1D_ptr(i)
      end do
    else if (rank == 2) then
      call mem_setptr(tgt2D_ptr, var_name, mem_path)
      call c_f_pointer(c_arr_ptr, src2D_ptr, shape(tgt2D_ptr))
      do j = 1, size(tgt2D_ptr, 2)
        do i = 1, size(tgt2D_ptr, 1)
          tgt2D_ptr(i, j) = src2D_ptr(i, j)
        end do
      end do
    else
      write (bmi_last_error, fmt_unsupported_rank) trim(var_name)
      call report_bmi_error(bmi_last_error)
      bmi_status = BMI_FAILURE
      return
    end if

    ! trigger event:
    call on_memory_set(var_name, mem_path, status)
    if (status /= 0) then
      ! something went terribly wrong here, aborting
      write (bmi_last_error, fmt_invalid_mem_access) trim(var_name)
      call report_bmi_error(bmi_last_error)
      bmi_status = BMI_FAILURE
      return
    end if

  end function set_value_int

  !> @brief Get the variable type as a string
  !!
  !! The type returned is that of a single element. Currently we
  !! support 'INTEGER' and 'DOUBLE'.  When the variable cannot
  !! be found, the string 'UNKNOWN' is assigned.
  !<
  function get_var_type(c_var_address, c_var_type) result(bmi_status) bind(C, name="get_var_type")
!DIR$ ATTRIBUTES DLLEXPORT :: get_var_type
    ! -- modules
    use ConstantsModule, only: LENMEMTYPE
    ! -- dummy variables
    character(kind=c_char), intent(in) :: c_var_address(*)             !< memory address string of the variable
    character(kind=c_char), intent(out) :: c_var_type(BMI_LENVARTYPE)  !< variable type as a string
    integer(kind=c_int) :: bmi_status                                   !< BMI status code
    ! -- local variables
    character(len=LENMEMPATH) :: mem_path
    character(len=LENVARNAME) :: var_name
    character(len=LENMEMTYPE) :: mem_type
    logical(LGP) :: valid

    bmi_status = BMI_SUCCESS

    call split_address(c_var_address, mem_path, var_name, valid)
    if (.not. valid) then
      bmi_status = BMI_FAILURE
      return
    end if

    call get_mem_type(var_name, mem_path, mem_type)
    c_var_type(1:len(trim(mem_type)) + 1) = string_to_char_array(trim(mem_type), len(trim(mem_type)))

    if (mem_type == 'UNKNOWN') then
      write (bmi_last_error, fmt_general_err) 'unknown memory type'
      call report_bmi_error(bmi_last_error)
      bmi_status = BMI_FAILURE
    end if

  end function get_var_type

  !> @brief Get the variable rank (non-BMI)
  !!
  !! In order to support multi-dimensional arrays, this function gives
  !! access to the rank of the array.
  !<
  function get_var_rank(c_var_address, c_var_rank) result(bmi_status) bind(C, name="get_var_rank")
    !DIR$ ATTRIBUTES DLLEXPORT :: get_var_rank
    ! -- dummy variables
    character(kind=c_char), intent(in) :: c_var_address(*) !< memory address string of the variable
    integer(kind=c_int), intent(out) :: c_var_rank          !< variable rank
    integer(kind=c_int) :: bmi_status                       !< BMI status code
    ! -- local variables
    character(len=LENMEMPATH) :: mem_path
    character(len=LENVARNAME) :: var_name
    logical(LGP) :: valid

    bmi_status = BMI_SUCCESS

    call split_address(c_var_address, mem_path, var_name, valid)
    if (.not. valid) then
      bmi_status = BMI_FAILURE
      return
    end if

    call get_mem_rank(var_name, mem_path, c_var_rank)
    if (c_var_rank == -1) then
      bmi_status = BMI_FAILURE
      return
    end if

  end function get_var_rank

  !> @brief Get the shape of the array for the variable (non-BMI)
  !!
  !! The shape is an integer array with size equal to the rank of
  !! the variable (see get_var_rank()) and values that give the
  !! length of the array in each dimension. The target shape array
  !! @p c_var_shape should be allocated before calling this routine.
  !!
  !! Note that the returned shape representation will has been converted
  !! to C-style.
  !<
  function get_var_shape(c_var_address, c_var_shape) result(bmi_status) bind(C, name="get_var_shape")
    !DIR$ ATTRIBUTES DLLEXPORT :: get_var_shape
    ! -- modules
    use ConstantsModule, only: MAXMEMRANK
    ! -- dummy variables
    character(kind=c_char), intent(in) :: c_var_address(*) !< memory address string of the variable
    integer(c_int), intent(inout) :: c_var_shape(*)         !< 1D array with the variable's shape
    integer(kind=c_int) :: bmi_status                       !< BMI status code
    ! -- local variables
    integer(I4B), dimension(MAXMEMRANK) :: var_shape
    integer(I4B) :: var_rank
    character(len=LENMEMPATH) :: mem_path
    character(len=LENVARNAME) :: var_name
    logical(LGP) :: valid

    bmi_status = BMI_SUCCESS

    call split_address(c_var_address, mem_path, var_name, valid)
    if (.not. valid) then
      bmi_status = BMI_FAILURE
      return
    end if

    var_shape = 0
    var_rank = 0
    call get_mem_rank(var_name, mem_path, var_rank)
    call get_mem_shape(var_name, mem_path, var_shape)
    if (var_shape(1) == -1 .or. var_rank == -1) then
      bmi_status = BMI_FAILURE
      return
    end if

    ! External calls to this BMI are assumed C style, so if the internal shape
    ! is (100,1) we get (100,1,undef) from the call get_mem_shape
    ! This we need to convert to C-style which should be (1,100).
    ! Hence, we reverse the array and drop undef:
    c_var_shape(1:var_rank) = var_shape(var_rank:1:-1)

  end function get_var_shape

end module mf6bmi


text:
      subroutine  sub_ilu2(neq,a,b,na,nb,ncon,nop
     *     ,irb,iirb,npvt,sorthm,dum1,dum2,piv)
!***********************************************************************
!  Copyright, 2004,  The  Regents  of the  University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
C**********************************************************************
CD1
CD1 PURPOSE
CD1
CD1 To perform incomplete lu factorization for 2 degree of freedom
CD1 problem.
CD1
C**********************************************************************
CD2
CD2 REVISION HISTORY
CD2
CD2 Revision                    ECD
CD2 Date         Programmer     Number  Comments
CD2
CD2 4-6-94       G. Zyvoloski   97      Initial implementation
CD2
CD2 $Log:   /pvcs.config/fehm90/src/sub_ilu2.f_a  $
!D2
!D2    Rev 2.5   06 Jan 2004 10:44:12   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 09:19:52   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!D2    Rev 2.3   14 Nov 2001 13:15:40   pvcs
!D2 FEHM Version 2.12, STN 10086-2.12-00
!D2
!D2    Rev 2.2   06 Jun 2001 13:27:48   pvcs
!D2 FEHM Version 2.11, STN 10086-2.11-00
!D2
!D2    Rev 2.1   30 Nov 2000 12:12:04   pvcs
!D2 FEHM Version 2.10, STN 10086-2.10-00
!D2
!D2    Rev 2.0   Fri May 07 14:46:56 1999   pvcs
!D2 FEHM Version 2.0, SC-194 (Fortran 90)
CD2
CD2    Rev 1.4   Thu Sep 12 08:26:40 1996   robinson
CD2 Prolog Changes
CD2
CD2    Rev 1.3   Fri Feb 02 12:41:26 1996   hend
CD2 Added Requirements Traceability
CD2
CD2    Rev 1.2   05/11/94 16:19:18   llt
CD2 bug fixes - gaz
CD2
CD2    Rev 1.1   03/18/94 16:05:42   gaz
CD2 Added solve_new and cleaned up memory management.
CD2
CD2    Rev 1.0   03/02/94 08:46:56   pvcs
CD2 original version
CD2
C**********************************************************************
CD3
CD3 INTERFACES
CD3
CD3 Formal Calling Parameters
CD3
CD3 Identifier   Type    Use     Description
CD3
CD3 neq          integer  I      Number of entries in the array for
CD3                                  each degree of freedom
CD3 a            real*8   I      Solution matrix
CD3 b            real*8   O      lu factorization matrix
CD3 na           integer  I      Pointer array in a matrix
CD3 nb           integer  I      Pointer array in b matrix
CD3 ncon         integer  I      Connectivity matrix for solution matrix
CD3 nop          integer  I      Connectivity matrix for factorization
CD3                                 matrix
CD3 irb          integer  I      Renumber array - inew=irb(iold)
CD3 iirb         integer  I      Renumber array - iold=iirb(inew)
CD3 npvt         integer  I      Pivot positions in nop
CD3 sorthm       real*8   I      Scratch storage array
CD3 dum1         real*8   I      Scratch storage array
CD3 dum2         real*8   I      Scratch storage array
CD3 piv          real*8   O      Array of pivots
CD3
CD3 Interface Tables
CD3
CD3 None
CD3
CD3 Files
CD3
CD3 NONE
CD3
C**********************************************************************
CD4
CD4 GLOBAL OBJECTS
CD4
CD4 Global Constants
CD4
CD4 NONE
CD4
CD4 Global Types
CD4
CD4 NONE
CD4
CD4 Global Variables
CD4
CD4 NONE
CD4
CD4 Global Subprograms
CD4
CD4 NONE
CD4
C**********************************************************************
CD5
CD5 LOCAL IDENTIFIERS
CD5
CD5 Local Constants
CD5
CD5 NONE
CD5
CD5 Local Types
CD5
CD5 NONE
CD5
CD5 Local variables
CD5
CD5 Identifier   Type        Description
CD5
CD5 neqm1        integer     neq-1
CD5 neqp1        integer     neq+1
CD5 na1          integer     Pointer value in a array
CD5 na2          integer     Pointer value in a array
CD5 na3          integer     Pointer value in a array
CD5 na4          integer     Pointer value in a array
CD5 nb1          integer     Pointer value in b array
CD5 nb2          integer     Pointer value in b array
CD5 nb3          integer     Pointer value in b array
CD5 nb4          integer     Pointer value in b array
CD5 i            integer     Do loop index
CD5 i1           integer     Do loop index
CD5 i2           integer     Do loop index
CD5 i3           integer     Do loop index
CD5 i4           integer     Do loop index
CD5 i5           integer     Do loop index
CD5 i6           integer     Do loop index
CD5 ip           integer     Pointer to reorder array
CD5 ik           integer     Do loop index
CD5 ikb          integer     Integer index
CD5 k            integer     Do loop index
CD5 ka           integer     Pointer in connectivity array
CD5 kb           integer     Pointer to reorder array
CD5 l            integer     Do loop index
CD5 npiv         integer     Pivot position
CD5 k1           integer     k+1
CD5 npivk1       integer     Index in pivot array
CD5 nz           integer     Pointer to multiplier of previous row
CD5 kd           integer     Integer index
CD5 jj           integer     Do loop index
CD5 nj           integer     Position in lu pointer array
CD5 npivnj       integer     Index in pivot array
CD5 ii           integer     Do loop index
CD5 nk           integer     Position in factorization connectivity
CD5                             array
CD5 ijk          integer     Do loop index
CD5 kj           integer     Integer index
CD5 idum         integer     Integer index
CD5 pp           real*8      Scratch storage space
CD5 a1           real*8      Variable in function definition
CD5 a2           real*8      Variable in function definition
CD5 a3           real*8      Variable in function definition
CD5 a4           real*8      Variable in function definition
CD5 b1           real*8      Value in b array
CD5 b2           real*8      Value in b array
CD5 b3           real*8      Value in b array
CD5 b4           real*8      Value in b array
CD5 c1           real*8      Value used in calculation
CD5 c2           real*8      Value used in calculation
CD5 c3           real*8      Value used in calculation
CD5 c4           real*8      Value used in calculation
CD5 p1           real*8      Value in b array
CD5 p2           real*8      Value in b array
CD5 p3           real*8      Value in b array
CD5 p4           real*8      Value in b array
CD5
CD5 Local Subprograms
CD5
CD5 Identifier    Type      Description
CD5
CD5 alm           real*8    Linear combination function
CD5 ali           real*8    Linear combination function
CD5
C**********************************************************************
CD6
CD6 ASSUMPTIONS AND LIMITATIONS
CD6
CD6 N/A
CD6
C**********************************************************************
CD7
CD7 SPECIAL COMMENTS
CD7
CD7 N/A
CD7
C**********************************************************************
CD8
CD8 REQUIREMENTS TRACEABILITY
CD8
CD8 3.2 Solve Linear Equation Set
CD8    3.2.1 Perform Incomplete Factorization
CD8 3.3 Provide Multiple Degree of Freedom Option
CD8
C**********************************************************************
CDA
CDA REFERENCES
CDA
CDA See GZSOLVE SRS, MMS, and SDD for documentation.
CDA
C**********************************************************************
CPS
CPS PSEUDOCODE
CPS
CPS BEGIN sub_ilu2
CPS
CPS Define integer parameters
CPS
CPS FOR each equation
CPS
CPS   Determine index parameters
CPS
CPS   For each position in dummy storage arrays
CPS     Initialize to zero
CPS   ENDFOR
CPS
CPS   FOR each nonzero position in solution array
CPS     Copy solution matrix values to dummy storage arrays
CPS   ENDFOR
CPS
CPS   FOR each nonzero position in lu factorization  array
CPS     Copy dummy storage array values into lu factorization array
CPS   ENDFOR
CPS
CPS ENDFOR each equation
CPS
CPS FOR each equation except the last
CPS
CPS   FOR each value in b array for this equation
CPS     Store b value in scratch storage array
CPS   ENDFOR
CPS
CPS   Identify pivot position
CPS   Invert diagonal element to form pivot
CPS
CPS   Save pivot value in arrays
CPS
CPS   FOR each nonzero element in current equation
CPS     Compute values in b array
CPS   ENDFOR
CPS
CPS   FOR each nonzero element of lu matrix
CPS     Copy value into dummy storage array
CPS   ENDFOR
CPS
CPS   FOR each nonzero element in lower diagonal part of lu matrix
CPS
CPS     Save term
CPS     Identify column number of that term
CPS
CPS     FOR each column number
CPS       Add linear combination (place multiplier in array sorthm)...
CPS       ... of previous row with row number equal to this column...
CPS       ... number (place in array dum)
CPS     ENDFOR
CPS
CPS   ENDFOR each nonzero element in lower diagonal part of lu matrix
CPS
CPS   FOR each nonzero position in lu factorization matrix
CPS     Copy value in dum to lu factorization matrix
CPS   ENDFOR
CPS
CPS   FOR each nonzero value of lu factoriation matrix below diagonal
CPS     Copy value in sorthm to lu factorization matrix
CPS   ENDFOR
CPS
CPS ENDFOR each equation except the last
CPS
CPS Compute pivot for last equation
CPS
CPS END sub_ilu2
CPS
C**********************************************************************
c
c two degree of freedom solver
c
c     call ilu factorization at preconditioner stage

      implicit none

      integer neq
      integer ncon(*),nop(*),irb(*),iirb(*),npvt(*)
      real*8 a(*),b(*),sorthm(*),piv(neq,4)
      real*8 pp(4)
      real*8 dum1(neq,2),dum2(neq,2)
      integer na(*),nb(*)
      real*8 a1, a2, a3, a4, b1, b2, b3, b4, c1, c2, c3, c4
      real*8 p1, p2, p3, p4
      integer neqm1, neqp1, na1, na2, na3, na4, nb1, nb2, nb3, nb4
      integer i, i1, i2, ip, i3, i4, ik, ikb, k, ka, kb, l, npiv
      integer k1, npivk1, i5, nz, i6, kd, jj, nj, npivnj, ii, nk
      integer kj, ijk, idum
      real*8 alm, ali
c
c***     linear algebra     ***
c
      alm(a1,a2,b1,b2)=a1*b1+a2*b2
      ali(a1,a2,a3,a4,b1)=b1/(a1*a4-a2*a3)
c
c     define some parameters
c
      neqm1=neq-1
      neqp1=neq+1
      na1=na(1)
      na2=na(2)
      na3=na(3)
      na4=na(4)
      nb1=nb(1)
      nb2=nb(2)
      nb3=nb(3)
      nb4=nb(4)
c
      do i=1,neq
         i1=nop(i)+1
         i2=nop(i+1)
         ip=irb(i)
         i3=ncon(ip)+1
         i4=ncon(ip+1)
         do ik=i1,i2
            ikb=nop(ik)
            dum1(ikb,1)=0.0d00
            dum1(ikb,2)=0.0d00
            dum2(ikb,1)=0.0d00
            dum2(ikb,2)=0.0d00
         enddo
         do k=i3,i4
            ka=ncon(k)
            kb=iirb(ka)
            dum1(kb,1)=a(k-neqp1+na1)
            dum1(kb,2)=a(k-neqp1+na2)
            dum2(kb,1)=a(k-neqp1+na3)
            dum2(kb,2)=a(k-neqp1+na4)
         enddo
         do l=i1,i2
            kb=nop(l)
            b(l-neqp1+nb1)=dum1(kb,1)
            b(l-neqp1+nb2)=dum1(kb,2)
            b(l-neqp1+nb3)=dum2(kb,1)
            b(l-neqp1+nb4)=dum2(kb,2)
         enddo
      enddo
c
c     start factorization
c
      do  k=1,neqm1
         npiv=npvt(k)
         do i=1,4
            pp(i)=b(npiv-neqp1+nb(i))
         enddo
         piv(k,1)=ali(pp(1),pp(2),pp(3),pp(4),pp(4))
         piv(k,2)=-ali(pp(1),pp(2),pp(3),pp(4),pp(2))
         piv(k,3)=-ali(pp(1),pp(2),pp(3),pp(4),pp(3))
         piv(k,4)=ali(pp(1),pp(2),pp(3),pp(4),pp(1))
c     divide row k and rhs by pivot
         i2=nop(k+1)
         do  i=npiv,i2
            b1=b(i-neqp1+nb(1))
            b2=b(i-neqp1+nb(2))
            b3=b(i-neqp1+nb(3))
            b4=b(i-neqp1+nb(4))
            b(i-neqp1+nb(1))=alm(piv(k,1),piv(k,2),b1,b3)
            b(i-neqp1+nb(2))=alm(piv(k,1),piv(k,2),b2,b4)
            b(i-neqp1+nb(3))=alm(piv(k,3),piv(k,4),b1,b3)
            b(i-neqp1+nb(4))=alm(piv(k,3),piv(k,4),b2,b4)
         enddo
         k1=k+1
c     factor row k1
         npivk1=npvt(k1)-1
         i5=nop(k1)+1
         nz=0
         i6=nop(k1+1)
         do  i=i5,i6
            kd=nop(i)
            dum1(kd,1)=b(i-neqp1+nb(1))
            dum1(kd,2)=b(i-neqp1+nb(2))
            dum2(kd,1)=b(i-neqp1+nb(3))
            dum2(kd,2)=b(i-neqp1+nb(4))
         enddo
         do jj=i5,npivk1
            nj=nop(jj)
            nz=nz+1
            c1=dum1(nj,1)
            c2=dum1(nj,2)
            c3=dum2(nj,1)
            c4=dum2(nj,2)
            sorthm(nz)=c1
            sorthm(nz+neq)=c2
            sorthm(nz+neq+neq)=c3
            sorthm(nz+neq+neq+neq)=c4
c     add linear combo of row nj
            i4=nop(nj+1)
            npivnj=npvt(nj)
            do ii=npivnj,i4
               nk=nop(ii)
               p1=b(ii-neqp1+nb(1))
               p2=b(ii-neqp1+nb(2))
               p3=b(ii-neqp1+nb(3))
               p4=b(ii-neqp1+nb(4))
               dum1(nk,1)=dum1(nk,1)-alm(c1,c2,p1,p3)
               dum1(nk,2)=dum1(nk,2)-alm(c1,c2,p2,p4)
               dum2(nk,1)=dum2(nk,1)-alm(c3,c4,p1,p3)
               dum2(nk,2)=dum2(nk,2)-alm(c3,c4,p2,p4)
            enddo
         enddo
         do i=npivk1+1,i6
            kj=nop(i)
            b(i-neqp1+nb(1))=dum1(kj,1)
            b(i-neqp1+nb(2))=dum1(kj,2)
            b(i-neqp1+nb(3))=dum2(kj,1)
            b(i-neqp1+nb(4))=dum2(kj,2)
         enddo
         nz=0
         do  ijk=i5,npivk1
            nz=nz+1
            idum=ijk-neqp1
            b(idum+nb(1))=sorthm(nz)
            b(idum+nb(2))=sorthm(nz+neq)
            b(idum+nb(3))=sorthm(nz+neq+neq)
            b(idum+nb(4))=sorthm(nz+neq+neq+neq)
         enddo
      enddo
      idum=nop(neqp1)-neqp1
      b1=b(idum+nb(1))
      b2=b(idum+nb(2))
      b3=b(idum+nb(3))
      b4=b(idum+nb(4))
      piv(neq,1)=ali(b1,b2,b3,b4,b4)
      piv(neq,2)=-ali(b1,b2,b3,b4,b2)
      piv(neq,3)=-ali(b1,b2,b3,b4,b3)
      piv(neq,4)=ali(b1,b2,b3,b4,b1)
c
      return
      end

ext: f

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
************************************************************************
* SUBROUTINE PLISU              ALL SYSTEMS                   97/01/22
* PURPOSE :
* EASY TO USE SUBROUTINE FOR LARGE-SCALE UNCONSTRAINED MINIMIZATION.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  II  IPAR(7)  INTEGER PAREMETERS:
*      IPAR(1)  MAXIMUM NUMBER OF ITERATIONS.
*      IPAR(2)  MAXIMUM NUMBER OF FUNCTION EVALUATIONS.
*      IPAR(3)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIS.
*      IPAR(4)  ESTIMATION INDICATOR. IPAR(4)=0-MINIMUM IS NOT
*         ESTIMATED. IPAR(4)=1-MINIMUM IS ESTIMATED BY THE VALUE
*         RPAR(6).
*      IPAR(5)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIS.
*      IPAR(6)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIS.
*      IPAR(7)  MAXIMUM NUMBER OF VARIABLE METRIC UPDATES.
*  RI  RPAR(9)  REAL PARAMETERS:
*      RPAR(1)  MAXIMUM STEPSIZE.
*      RPAR(2)  TOLERANCE FOR THE CHANGE OF VARIABLES.
*      RPAR(3)  TOLERANCE FOR THE CHANGE OF FUNCTION VALUES.
*      RPAR(4)  TOLERANCE FOR THE FUNCTION FALUE.
*      RPAR(5)  TOLERANCE FOR THE GRADIENT NORM.
*      RPAR(6)  ESTIMATION OF THE MINIMUM FUNCTION VALUE.
*      RPAR(7)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIS.
*      RPAR(8)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIS.
*      RPAR(9)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIS.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE.
*  II  IPRNT  PRINT SPECIFICATION. IPRNT=0-NO PRINT.
*         ABS(IPRNT)=1-PRINT OF FINAL RESULTS.
*         ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS.
*         IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL
*         RESULTS.
*  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION.
*         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN
*                   MTESX (USUALLY TWO) SUBSEQUENT ITERATIONS.
*         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN
*                   MTESF (USUALLY TWO) SUBSEQUENT ITERATIONS.
*         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB.
*         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG.
*         ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED,
*                   BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE.
*         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV.
*         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED.
*
* VARIABLES IN COMMON /STAT/ (STATISTICS) :
*  IO  NRES  NUMBER OF RESTARTS.
*  IO  NDEC  NUMBER OF MATRIX DECOMPOSITIONS.
*  IO  NIN  NUMBER OF INNER ITERATIONS.
*  IO  NIT  NUMBER OF ITERATIONS.
*  IO  NFV  NUMBER OF FUNCTION EVALUATIONS.
*  IO  NFG  NUMBER OF GRADIENT EVALUATIONS.
*  IO  NFH  NUMBER OF HESSIAN EVALUATIONS.
*
* SUBPROGRAMS USED :
*  S   PLIS  LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG
*         RECURRENCES.
*
* EXTERNAL SUBROUTINES :
*  SE  OBJ  COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION.
*         CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER
*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS
*         THE VALUE OF THE OBJECTIVE FUNCTION.
*  SE  DOBJ  COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION.
*         CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER
*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF)
*         IS THE GRADIENT OF THE OBJECTIVE FUNCTION.
*
      SUBROUTINE PLISU(NF,X,IPAR,RPAR,F,GMAX,IPRNT,ITERM)
      INTEGER NF,IPAR(7),IPRNT,ITERM
      DOUBLE PRECISION X(*),RPAR(9),F,GMAX
      INTEGER MF,NB,LGF,LS,LXO,LGO,LUO,LVO
      INTEGER NRES,NDEC,NIN,NIT,NFV,NFG,NFH
      COMMON /STAT/ NRES,NDEC,NIN,NIT,NFV,NFG,NFH
      DOUBLE PRECISION RA(:)
      ALLOCATABLE RA
      MF=IPAR(7)
      IF (MF.LE.0) MF=10
      ALLOCATE (RA(2*NF+2*NF*MF+2*MF))
      NB=0
*
*     POINTERS FOR AUXILIARY ARRAYS
*
      LGF=1
      LS=LGF+NF
      LXO=LS+NF
      LGO=LXO+NF*MF
      LUO=LGO+NF*MF
      LVO=LUO+MF
      CALL PLIS(NF,NB,X,IPAR,RA,RA,RA(LGF),RA(LS),RA(LXO),RA(LGO),
     & RA(LUO),RA(LVO),RPAR(1),RPAR(2),RPAR(3),RPAR(4),RPAR(5),RPAR(6),
     & GMAX,F,IPAR(1),IPAR(2),IPAR(4),MF,IPRNT,ITERM)
      DEALLOCATE (RA)
      RETURN
      END
************************************************************************
* SUBROUTINE PLISS              ALL SYSTEMS                   97/01/22
* PURPOSE :
* EASY TO USE SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE
*         X(I) IS UNBOUNDED. IX(I)=1-LOWER BOUND XL(I).LE.X(I).
*         IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND
*         XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  II  IPAR(7)  INTEGER PAREMETERS:
*      IPAR(1)  MAXIMUM NUMBER OF ITERATIONS.
*      IPAR(2)  MAXIMUM NUMBER OF FUNCTION EVALUATIONS.
*      IPAR(3)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIS.
*      IPAR(4)  ESTIMATION INDICATOR. IPAR(4)=0-MINIMUM IS NOT
*         ESTIMATED. IPAR(4)=1-MINIMUM IS ESTIMATED BY THE VALUE
*         RPAR(6).
*      IPAR(5)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIS.
*      IPAR(6)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIS.
*      IPAR(7)  MAXIMUM NUMBER OF VARIABLE METRIC UPDATES.
*  RI  RPAR(9)  REAL PARAMETERS:
*      RPAR(1)  MAXIMUM STEPSIZE.
*      RPAR(2)  TOLERANCE FOR THE CHANGE OF VARIABLES.
*      RPAR(3)  TOLERANCE FOR THE CHANGE OF FUNCTION VALUES.
*      RPAR(4)  TOLERANCE FOR THE FUNCTION FALUE.
*      RPAR(5)  TOLERANCE FOR THE GRADIENT NORM.
*      RPAR(6)  ESTIMATION OF THE MINIMUM FUNCTION VALUE.
*      RPAR(7)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIS.
*      RPAR(8)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIS.
*      RPAR(9)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIS.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE.
*  II  IPRNT  PRINT SPECIFICATION. IPRNT=0-NO PRINT.
*         ABS(IPRNT)=1-PRINT OF FINAL RESULTS.
*         ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS.
*         IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL
*         RESULTS.
*  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION.
*         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN
*                   MTESX (USUALLY TWO) SUBSEQUENT ITERATIONS.
*         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN
*                   MTESF (USUALLY TWO) SUBSEQUENT ITERATIONS.
*         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB.
*         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG.
*         ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED,
*                   BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE.
*         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV.
*         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED.
*
* VARIABLES IN COMMON /STAT/ (STATISTICS) :
*  IO  NRES  NUMBER OF RESTARTS.
*  IO  NDEC  NUMBER OF MATRIX DECOMPOSITIONS.
*  IO  NIN  NUMBER OF INNER ITERATIONS.
*  IO  NIT  NUMBER OF ITERATIONS.
*  IO  NFV  NUMBER OF FUNCTION EVALUATIONS.
*  IO  NFG  NUMBER OF GRADIENT EVALUATIONS.
*  IO  NFH  NUMBER OF HESSIAN EVALUATIONS.
*
* SUBPROGRAMS USED :
*  S   PLIS  LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG
*         RECURRENCES.
*
* EXTERNAL SUBROUTINES :
*  SE  OBJ  COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION.
*         CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER
*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS
*         THE VALUE OF THE OBJECTIVE FUNCTION.
*  SE  DOBJ  COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION.
*         CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER
*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF)
*         IS THE GRADIENT OF THE OBJECTIVE FUNCTION.
*
      SUBROUTINE PLISS(NF,X,IX,XL,XU,IPAR,RPAR,F,GMAX,IPRNT,ITERM)
      INTEGER NF,IX(*),IPAR(7),IPRNT,ITERM
      DOUBLE PRECISION X(*),XL(*),XU(*),RPAR(9),F,GMAX
      INTEGER MF,NB,LGF,LS,LXO,LGO,LUO,LVO
      INTEGER NRES,NDEC,NIN,NIT,NFV,NFG,NFH
      COMMON /STAT/ NRES,NDEC,NIN,NIT,NFV,NFG,NFH
      DOUBLE PRECISION RA(:)
      ALLOCATABLE RA
      MF=IPAR(7)
      IF (MF.LE.0) MF=10
      ALLOCATE (RA(2*NF+2*NF*MF+2*MF))
      NB=1
*
*     POINTERS FOR AUXILIARY ARRAYS
*
      LGF=1
      LS=LGF+NF
      LXO=LS+NF
      LGO=LXO+NF*MF
      LUO=LGO+NF*MF
      LVO=LUO+MF
      CALL PLIS(NF,NB,X,IX,XL,XU,RA(LGF),RA(LS),RA(LXO),RA(LGO),
     & RA(LUO),RA(LVO),RPAR(1),RPAR(2),RPAR(3),RPAR(4),RPAR(5),RPAR(6),
     & GMAX,F,IPAR(1),IPAR(2),IPAR(4),MF,IPRNT,ITERM)
      DEALLOCATE (RA)
      RETURN
      END
************************************************************************
* SUBROUTINE PLIS               ALL SYSTEMS                   01/09/22
* PURPOSE :
* GENERAL SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION THAT
* USE THE LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG
* RECURRENCES.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  NB  CHOICE OF SIMPLE BOUNDS. NB=0-SIMPLE BOUNDS SUPPRESSED.
*         NB>0-SIMPLE BOUNDS ACCEPTED.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE
*         X(I) IS UNBOUNDED. IX(I)=1-LOVER BOUND XL(I).LE.X(I).
*         IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND
*         XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RO  GF(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RO  S(NF)  DIRECTION VECTOR.
*  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE.
*  RI  GO(NF)  GRADIENTS DIFFERENCE.
*  RA  UO(NF)  AUXILIARY VECTOR.
*  RA  VO(NF)  AUXILIARY VECTOR.
*  RI  XMAX  MAXIMUM STEPSIZE.
*  RI  TOLX  TOLERANCE FOR CHANGE OF VARIABLES.
*  RI  TOLF  TOLERANCE FOR CHANGE OF FUNCTION VALUES.
*  RI  TOLB  TOLERANCE FOR THE FUNCTION VALUE.
*  RI  TOLG  TOLERANCE FOR THE GRADIENT NORM.
*  RI  FMIN  ESTIMATION OF THE MINIMUM FUNCTION VALUE.
*  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  II  MIT  MAXIMUM NUMBER OF ITERATIONS.
*  II  MFV  MAXIMUM NUMBER OF FUNCTION EVALUATIONS.
*  II  IEST  ESTIMATION INDICATOR. IEST=0-MINIMUM IS NOT ESTIMATED.
*         IEST=1-MINIMUM IS ESTIMATED BY THE VALUE FMIN.
*  II  MF  NUMBER OF LIMITED MEMORY STEPS.
*  II  IPRNT  PRINT SPECIFICATION. IPRNT=0-NO PRINT.
*         ABS(IPRNT)=1-PRINT OF FINAL RESULTS.
*         ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS.
*         IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL
*         RESULTS.
*  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION.
*         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN
*                   MTESX (USUALLY TWO) SUBSEQUEBT ITERATIONS.
*         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN
*                   MTESF (USUALLY TWO) SUBSEQUEBT ITERATIONS.
*         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB.
*         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG.
*         ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED,
*                   BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE.
*         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV.
*         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED.
*
* VARIABLES IN COMMON /STAT/ (STATISTICS) :
*  IO  NRES  NUMBER OF RESTARTS.
*  IO  NDEC  NUMBER OF MATRIX DECOMPOSITION.
*  IO  NIN  NUMBER OF INNER ITERATIONS.
*  IO  NIT  NUMBER OF ITERATIONS.
*  IO  NFV  NUMBER OF FUNCTION EVALUATIONS.
*  IO  NFG  NUMBER OF GRADIENT EVALUATIONS.
*  IO  NFH  NUMBER OF HESSIAN EVALUATIONS.
*
* SUBPROGRAMS USED :
*  S   PCBS04  ELIMINATION OF BOX CONSTRAINT VIOLATIONS.
*  S   PS1L01  STEPSIZE SELECTION USING LINE SEARCH.
*  S   PYADC0  ADDITION OF A BOX CONSTRAINT.
*  S   PYFUT1  TEST ON TERMINATION.
*  S   PYRMC0  DELETION OF A BOX CONSTRAINT.
*  S   PYTRCD  COMPUTATION OF PROJECTED DIFFERENCES FOR THE VARIABLE METRIC
*         UPDATE.
*  S   PYTRCG  COMPUTATION OF THE PROJECTED GRADIENT.
*  S   PYTRCS  COMPUTATION OF THE PROJECTED DIRECTION VECTOR.
*  S   MXDRCB BACKWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION
*         OF THE VECTOR X BY AN IMPLICIT BFGS UPDATE.
*  S   MXDRCF FORWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION
*         OF THE VECTOR X BY AN IMPLICIT BFGS UPDATE.
*  S   MXDRSU SHIFT OF COLUMNS OF THE RECTANGULAR MATRICES A AND B.
*         SHIFT OF ELEMENTS OF THE VECTOR U. THESE SHIFTS ARE USED IN
*         THE LIMITED MEMORY BFGS METHOD.
*  S   MXUDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  RF  MXUDOT  DOT PRODUCT OF TWO VECTORS.
*  S   MXUNEG  COPYING OF A VECTOR WITH CHANGE OF THE SIGN.
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVSCL  SCALING OF A VECTOR.
*
* EXTERNAL SUBROUTINES :
*  SE  OBJ  COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION.
*         CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER
*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS
*         THE VALUE OF THE OBJECTIVE FUNCTION.
*  SE  DOBJ  COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION.
*         CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER
*         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF)
*         IS THE GRADIENT OF THE OBJECTIVE FUNCTION.
*
* METHOD :
* LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG
* RECURRENCES.
*
      SUBROUTINE PLIS(NF,NB,X,IX,XL,XU,GF,S,XO,GO,UO,VO,XMAX,TOLX,
     & TOLF,TOLB,TOLG,FMIN,GMAX,F,MIT,MFV,IEST,MF,IPRNT,ITERM)
      INTEGER NF,NB,IX(*),MIT,MFV,IEST,MF,IPRNT,ITERM
      DOUBLE PRECISION X(*),XL(*),XU(*),GF(*),S(*),XO(*),GO(*),UO(*),
     & VO(*),TOLX,TOLF,TOLG,TOLB,FMIN,XMAX,GMAX,F
      INTEGER ITERD,ITERS,KD,LD,NTESX,NTESF,MTESX,MTESF,MRED,KIT,
     & IREST,KBF,MES,MES1,MES2,MES3,MAXST,ISYS,ITES,INITS,KTERS,
     & IRES1,IRES2,INEW,IOLD,I,N,MFG,K,NRED
      DOUBLE PRECISION R,RO,RP,FO,FP,P,PO,PP,GNORM,SNORM,RMIN,RMAX,
     & UMAX,FMAX,DMAX,ETA0,ETA9,EPS8,EPS9,ALF1,ALF2,PAR1,PAR2,A,B,
     & TOLD,TOLS,TOLP
      DOUBLE PRECISION MXUDOT
      INTEGER NRES,NDEC,NIN,NIT,NFV,NFG,NFH
      COMMON /STAT/ NRES,NDEC,NIN,NIT,NFV,NFG,NFH
      IF (ABS(IPRNT).GT.1) WRITE(6,'(1X,''ENTRY TO PLIS :'')')
*
*     INITIATION
*
      KBF=0
      IF (NB.GT.0) KBF=2
      NRES=0
      NDEC=0
      NIN=0
      NIT=0
      NFV=0
      NFG=0
      NFH=0
      ISYS=0
      ITES=1
      MTESX=2
      MTESF=2
      INITS=2
      ITERM=0
      ITERD=0
      ITERS=2
      KTERS=3
      IREST=0
      IRES1=999
      IRES2=0
      MRED=10
      MES=4
      MES1=2
      MES2=2
      MES3=2
      ETA0=1.0D-15
      ETA9=1.0D 120
      EPS8=1.00D 0
      EPS9=1.00D-8
      ALF1=1.0D-10
      ALF2=1.0D 10
      RMAX=ETA9
      DMAX=ETA9
      FMAX=1.0D 20
      IF (IEST.LE.0) FMIN=-1.0D 60
      IF (IEST.GT.0) IEST=1
      IF (XMAX.LE.0.0D 0) XMAX=1.0D 16
      IF (TOLX.LE.0.0D 0) TOLX=1.0D-16
      IF (TOLF.LE.0.0D 0) TOLF=1.0D-14
      IF (TOLG.LE.0.0D 0) TOLG=1.0D-6
      IF (TOLB.LE.0.0D 0) TOLB=FMIN+1.0D-16
      TOLD=1.0D-4
      TOLS=1.0D-4
      TOLP=0.8D 0
      IF (MIT.LE.0) MIT=9000
      IF (MFV.LE.0) MFV=9000
      MFG=MFV
      KD= 1
      LD=-1
      KIT=-(IRES1*NF+IRES2)
      FO=FMIN
*
*     INITIAL OPERATIONS WITH SIMPLE BOUNDS
*
      IF (KBF.GT.0) THEN
      DO 2 I = 1,NF
      IF ((IX(I).EQ.3.OR.IX(I).EQ.4) .AND. XU(I).LE.XL(I)) THEN
      XU(I) = XL(I)
      IX(I) = 5
      ELSE IF (IX(I).EQ.5 .OR. IX(I).EQ.6) THEN
      XL(I) = X(I)
      XU(I) = X(I)
      IX(I) = 5
      END IF
    2 CONTINUE
      CALL PCBS04(NF,X,IX,XL,XU,EPS9,KBF)
      CALL PYADC0(NF,N,X,IX,XL,XU,INEW)
      END IF
      IF (ITERM.NE.0) GO TO 11190
      CALL OBJ(NF,X,F)
      NFV=NFV+1
      CALL DOBJ(NF,X,GF)
      NFG=NFG+1
11120 CONTINUE
      CALL PYTRCG(NF,NF,IX,GF,UMAX,GMAX,KBF,IOLD)
      IF (ABS(IPRNT).GT.1)
     & WRITE (6,'(1X,''NIT='',I5,2X,''NFV='',I5,2X,''NFG='',I5,2X,
     & ''F='', G16.9,2X,''G='',E10.3)') NIT,NFV,NFG,F,GMAX
      CALL PYFUT1(NF,F,FO,UMAX,GMAX,DMAX,TOLX,TOLF,TOLB,TOLG,KD,
     & NIT,KIT,MIT,NFV,MFV,NFG,MFG,NTESX,MTESX,NTESF,MTESF,ITES,
     & IRES1,IRES2,IREST,ITERS,ITERM)
      IF (ITERM.NE.0) GO TO 11190
      IF (KBF.GT.0.AND.RMAX.GT.0.0D 0) THEN
      CALL PYRMC0(NF,N,IX,GF,EPS8,UMAX,GMAX,RMAX,IOLD,IREST)
      END IF
11130 CONTINUE
*
*     DIRECTION DETERMINATION
*
      GNORM=SQRT(MXUDOT(NF,GF,GF,IX,KBF))
      IF (IREST.NE.0) GO TO 12620
      K=MIN(NIT-KIT,MF)
      IF (K.LE.0) THEN
      IREST=MAX(IREST,1)
      GO TO 12620
      END IF
*
*     DETERMINATION OF THE PARAMETER B
*
      B=MXUDOT(NF,XO,GO,IX,KBF)
      IF (B.LE.0.0D 0) THEN
      IREST=MAX(IREST,1)
      GO TO 12620
      END IF
      UO(1)=1.0D 0/B
      CALL MXUNEG(NF,GF,S,IX,KBF)
      CALL MXDRCB(NF,K,XO,GO,UO,VO,S,IX,KBF)
      A=MXUDOT(NF,GO,GO,IX,KBF)
      IF (A.GT.0.0D 0) THEN
      CALL MXVSCL(NF,B/A,S,S)
      END IF
      CALL MXDRCF(NF,K,XO,GO,UO,VO,S,IX,KBF)
      SNORM=SQRT(MXUDOT(NF,S,S,IX,KBF))
      K=MIN(K+1,MF)
      CALL MXDRSU(NF,K,XO,GO,UO)
12620 CONTINUE
      ITERD=0
      IF (IREST.NE.0) THEN
*
*     STEEPEST DESCENT DIRECTION
*
      CALL MXUNEG(NF,GF,S,IX,KBF)
      SNORM=GNORM
      IF (KIT.LT.NIT) THEN
        NRES=NRES+1
        KIT = NIT
      ELSE
        ITERM=-10
        IF (ITERS.LT.0) ITERM=ITERS-5
      END IF
      END IF
*
*     TEST ON DESCENT DIRECTION AND PREPARATION OF LINE SEARCH
*
      IF (KD.GT.0) P=MXUDOT(NF,GF,S,IX,KBF)
      IF (ITERD.LT.0) THEN
        ITERM=ITERD
      ELSE
*
*     TEST ON DESCENT DIRECTION
*
      IF (SNORM.LE.0.0D 0) THEN
        IREST=MAX(IREST,1)
      ELSE IF (P+TOLD*GNORM*SNORM.LE.0.0D 0) THEN
        IREST=0
      ELSE
*
*     UNIFORM DESCENT CRITERION
*
      IREST=MAX(IREST,1)
      END IF
      IF (IREST.EQ.0) THEN
*
*     PREPARATION OF LINE SEARCH
*
        NRED = 0
        RMIN=ALF1*GNORM/SNORM
        RMAX=MIN(ALF2*GNORM/SNORM,XMAX/SNORM)
      END IF
      END IF
      IF (ITERM.NE.0) GO TO 11190
      IF (IREST.NE.0) GO TO 11130
      CALL PYTRCS(NF,X,IX,XO,XL,XU,GF,GO,S,RO,FP,FO,F,PO,P,RMAX,ETA9,
     & KBF)
      IF (RMAX.EQ.0.0D 0) GO TO 11175
11170 CONTINUE
      CALL PS1L01(R,RP,F,FO,FP,P,PO,PP,FMIN,FMAX,RMIN,RMAX,
     & TOLS,TOLP,PAR1,PAR2,KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,
     & INITS,ITERS,KTERS,MES,ISYS)
      IF (ISYS.EQ.0) GO TO 11174
      CALL MXUDIR(NF,R,S,XO,X,IX,KBF)
      CALL PCBS04(NF,X,IX,XL,XU,EPS9,KBF)
      CALL OBJ(NF,X,F)
      NFV=NFV+1
      CALL DOBJ(NF,X,GF)
      NFG=NFG+1
      P=MXUDOT(NF,GF,S,IX,KBF)
      GO TO 11170
11174 CONTINUE
      IF (ITERS.LE.0) THEN
      R=0.0D 0
      F=FO
      P=PO
      CALL MXVCOP(NF,XO,X)
      CALL MXVCOP(NF,GO,GF)
      IREST=MAX(IREST,1)
      LD=KD
      GO TO 11130
      END IF
      CALL PYTRCD(NF,X,IX,XO,GF,GO,R,F,FO,P,PO,DMAX,KBF,KD,LD,ITERS)
11175 CONTINUE
      IF (KBF.GT.0) THEN
      CALL MXVINE(NF,IX)
      CALL PYADC0(NF,N,X,IX,XL,XU,INEW)
      END IF
      GO TO 11120
11190 CONTINUE
      IF (IPRNT.GT.1.OR.IPRNT.LT.0)
     & WRITE(6,'(1X,''EXIT FROM PLIS :'')')
      IF (IPRNT.NE.0)
     & WRITE (6,'(1X,''NIT='',I5,2X,''NFV='',I5,2X,''NFG='',I5,2X,
     & ''F='', G16.9,2X,''G='',E10.3,2X,''ITERM='',I3)') NIT,NFV,NFG,
     & F,GMAX,ITERM
      IF (IPRNT.LT.0)
     & WRITE (6,'(1X,''X='',5(G14.7,1X):/(3X,5(G14.7,1X)))')
     & (X(I),I=1,NF)
      RETURN
      END

ext: for

max_stars_repo_name: DenisSch/svca

max_issues_repo_name: DenisSch/svca

max_forks_repo_name: DenisSch/svca


text:
        SUBROUTINE ENDGET( BLOCK )
      INCLUDE 'DSIOF.COM'
      INCLUDE 'XNSTRN.COM'
        INTEGER BLOCK( 15 )
        NAME = BLOCK( 1 )
        CALL DSGEFL
        NWORDS = BLOCK( 11 )
        NELM = IAND( IBASE( INDCBP-2 ), MASKH2 )
        INDCBP = INDCBP + NELM*NWORDS + BLOCK(3)*2
        CALL DSSDCB
        RETURN
        END


!
!                Parallel Sparse BLAS  version 3.5
!      (C) Copyright 2006-2018
!        Salvatore Filippone
!        Alfredo Buttari
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!
!
! File: psb_dspsm.f90
!
!
! Subroutine: psb_dsps_vect
!  Performs one of the distributed matrix-vector operations
!
!     Y := alpha * Pr * A-1 * Pc * X + beta * Y,   or
!
!     Y := alpha * D * Pr * A-1 * Pc * X + beta * Y, or
!
!     Y := alpha * Pr * A-1 * Pc * D * X + beta * Y,  or
!
!     Y := alpha * Pr * A-T * Pc * X + beta * Y,   or
!
!     Y := alpha * D * Pr * A-T * Pc * X + beta * Y, or
!
!     Y := alpha * Pr * A-T * Pc * D * X + beta * Y,  or
!
!  X is a distributed
!  vector and T is a M-by-M distributed triangular matrix.
!
!
! Arguments:
!    alpha   -  real.               The scalar alpha.
!    a       -  type(psb_dspmat_type). The sparse matrix containing A.
!    x       -  type(psb_d_vect_type) The input vector containing the entries of ( X ).
!    beta    -  real                The scalar beta.
!    y       -  type(psb_d_vect_type) The input vector containing the entries of ( Y ).
!    desc_a  -  type(psb_desc_type).   The communication descriptor.
!    info    -  integer.               Return code
!    trans   -  character(optional).   Whether A or A'. If not present 'N' is assumed.
!    scale   -  character(optional).   Specify some type of operation with
!                                      the diagonal matrix D.
!    choice  -  integer(optional).     The kind of update to perform on overlap elements.
!    d(:)    -  real, optional      Matrix for diagonal scaling.
!    work(:) -  real, optional      Working area.
!
subroutine  psb_dspsv_vect(alpha,a,x,beta,y,desc_a,info,&
     & trans, scale, choice, diag, work)
  use psb_base_mod, psb_protect_name => psb_dspsv_vect
  use psi_mod
  implicit none

  real(psb_dpk_), intent(in)           :: alpha, beta
  type(psb_d_vect_type), intent(inout)    :: x
  type(psb_d_vect_type), intent(inout)    :: y
  type(psb_dspmat_type), intent(inout)    :: a
  type(psb_desc_type), intent(in)         :: desc_a
  integer(psb_ipk_), intent(out)                    :: info
  type(psb_d_vect_type), intent(inout), optional  :: diag
  real(psb_dpk_), optional, target, intent(inout) :: work(:)
  character, intent(in), optional         :: trans, scale
  integer(psb_ipk_), intent(in), optional           :: choice

  ! locals
  type(psb_ctxt_type) :: ctxt
  integer(psb_ipk_) :: np, me, &
       & err_act, iix, jjx, ia, ja, iia, jja, lldx,lldy, choice_,&
       & ix, iy, ik, jx, jy, i, lld,&
       & m, nrow, ncol, liwork, llwork, iiy, jjy, idx, ndm

  character                :: lscale
  integer(psb_ipk_), parameter       :: nb=4
  real(psb_dpk_),pointer :: iwork(:), xp(:), yp(:)
  character                :: itrans
  character(len=20)        :: name, ch_err
  logical                  :: aliw

  name='psb_dspsv_vect'
  info=psb_success_
  call psb_erractionsave(err_act)
  if  (psb_errstatus_fatal()) then
    info = psb_err_internal_error_ ;    goto 9999
  end if

  ctxt=desc_a%get_context()

  call psb_info(ctxt, me, np)
  if (np == -1) then
    info = psb_err_context_error_
    call psb_errpush(info,name)
    goto 9999
  endif

  if (.not.allocated(x%v)) then
    info = psb_err_invalid_vect_state_
    call psb_errpush(info,name)
    goto 9999
  endif
  if (.not.allocated(y%v)) then
    info = psb_err_invalid_vect_state_
    call psb_errpush(info,name)
    goto 9999
  endif

  if (present(choice)) then
    choice_ = choice
  else
    choice_ = psb_avg_
  endif

  if (present(scale)) then
    lscale = psb_toupper(scale)
  else
    lscale = 'U'
  endif

  if (present(trans)) then
    itrans = psb_toupper(trans)
    if((itrans == 'N').or.(itrans == 'T').or.(itrans == 'C')) then
      ! Ok
    else
      info = psb_err_iarg_invalid_value_
      call psb_errpush(info,name)
      goto 9999
    end if
  else
    itrans = 'N'
  endif

  m    = desc_a%get_global_rows()
  nrow = desc_a%get_local_rows()
  ncol = desc_a%get_local_cols()
  lldx = x%get_nrows()
  lldy = y%get_nrows()
  if ((info == 0).and.(lldx<ncol)) call x%reall(ncol,info)
  if ((info == 0).and.(lldy<ncol)) call y%reall(ncol,info)

  if (psb_errstatus_fatal()) then
    info=psb_err_from_subroutine_
    ch_err='reall'
    call psb_errpush(info,name,a_err=ch_err)
    goto 9999
  end if

  iwork => null()
  ! check for presence/size of a work area
  liwork= 2*ncol

  if (present(work)) then
    if (size(work) >= liwork) then
      aliw =.false.
    else
      aliw=.true.
    endif
  else
    aliw=.true.
  end if

  if (aliw) then
    allocate(iwork(liwork),stat=info)
    if(info /= psb_success_) then
      info=psb_err_from_subroutine_
      ch_err='psb_realloc'
      call psb_errpush(info,name,a_err=ch_err)
      goto 9999
    end if
  else
    iwork => work
  endif

  iwork(1)=0.d0

  ! Perform local triangular system solve
  if (present(diag)) then
    call a%spsm(alpha,x,beta,y,info,scale=scale,d=diag,trans=trans)
  else
    call a%spsm(alpha,x,beta,y,info,scale=scale,trans=trans)
  end if
  if(info /= psb_success_) then
    info = psb_err_from_subroutine_
    ch_err='dcssm'
    call psb_errpush(info,name,a_err=ch_err)
    goto 9999
  end if

  ! update overlap elements
  if (choice_ > 0) then
    call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),&
         & done,y%v,desc_a,iwork,info,data=psb_comm_ovr_)

    if (info == psb_success_) call psi_ovrl_upd(y%v,desc_a,choice_,info)
    if (info /= psb_success_) then
      call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner updates')
      goto 9999
    end if
  end if

  if (aliw) deallocate(iwork)

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(ctxt,err_act)

  return
end subroutine psb_dspsv_vect
!
! Subroutine: psb_dspsm
!  Performs one of the distributed matrix-vector operations
!
!     sub( Y ) := alpha * Pr * A-1 * Pc *sub( X ) + beta * sub (Y ),   or
!
!     sub( Y ) := alpha * D * Pr * A-1 * Pc * sub( X ) + beta * sub (Y ), or
!
!     sub( Y ) := alpha * Pr * A-1 * Pc * D * sub( X ) + beta * sub (Y ),  or
!
!     sub( Y ) := alpha * Pr * A-T * Pc * sub( X ) + beta * sub (Y ),   or
!
!     sub( Y ) := alpha * D * Pr * A-T * Pc * sub( X ) + beta * sub (Y ), or
!
!     sub( Y ) := alpha * Pr * A-T * Pc * D * sub( X ) + beta * sub (Y ),  or
!
!  where :
!
!        sub( X ) denotes X(1:M,JX:JX+K-1),
!
!        sub( Y ) denotes Y(1:M,JY:JY+K-1).
!
!  sub( X ) is a distributed
!  vector and T is a M-by-M distributed triangular matrix.
!
! Arguments:
!    alpha   -  real.               The scalar alpha.
!    a       -  type(psb_dspmat_type). The sparse matrix containing A.
!    x(:,:)  -  real                The input vector containing the entries of ( X ).
!    beta    -  real                The scalar beta.
!    y(:,:)  -  real                The input vector containing the entries of ( Y ).
!    desc_a  -  type(psb_desc_type).   The communication descriptor.
!    info    -  integer.               Return code
!    trans   -  character(optional).   Whether A or A'. If not present 'N' is assumed.
!    scale   -  character(optional).   Specify some type of operation with
!                                      the diagonal matrix D.
!    choice  -  integer(optional).     The kind of update to perform on overlap elements.
!    d(:)    -  real, optional      Matrix for diagonal scaling.
!    k       -  integer(optional).     The number of right-hand sides.
!    jx      -  integer(optional).     The column offset for ( X ). Default: 1
!    jy      -  integer(optional).     The column offset for ( Y ). Default: 1
!    work(:) -  real, optional      Working area.
!
subroutine  psb_dspsm(alpha,a,x,beta,y,desc_a,info,&
     & trans, scale, choice, diag, k, jx, jy, work)
  use psb_base_mod, psb_protect_name => psb_dspsm
  use psi_mod
  implicit none

  real(psb_dpk_), intent(in)              :: alpha, beta
  real(psb_dpk_), intent(in), target      :: x(:,:)
  real(psb_dpk_), intent(inout), target   :: y(:,:)
  type (psb_dspmat_type), intent(in)        :: a
  type(psb_desc_type), intent(in)           :: desc_a
  integer(psb_ipk_), intent(out)                      :: info
  real(psb_dpk_), intent(in), optional, target      :: diag(:)
  real(psb_dpk_), optional, target, intent(inout)   :: work(:)
  character, intent(in), optional           :: trans, scale
  integer(psb_ipk_), intent(in), optional             :: choice
  integer(psb_ipk_), intent(in), optional             :: k, jx, jy

  ! locals
  type(psb_ctxt_type) :: ctxt
  integer(psb_ipk_) :: np, me,&
       & err_act, iix, jjx, iia, jja, lldx,lldy, choice_,&
       & ik, i, lld, nrow, ncol, liwork, llwork, iiy, jjy, idx, ndm

  integer(psb_lpk_) :: ix, ijx, iy, ijy, m, n, ia, ja, lik
  character                :: lscale
  integer(psb_ipk_), parameter  :: nb=4
  real(psb_dpk_),pointer :: iwork(:), xp(:,:), yp(:,:), id(:)
  character                :: itrans
  character(len=20)        :: name, ch_err
  logical                  :: aliw

  name='psb_dspsm'
  info=psb_success_
  call psb_erractionsave(err_act)
  if (psb_errstatus_fatal()) then
    info = psb_err_internal_error_ ;    goto 9999
  end if

  ctxt=desc_a%get_context()

  call psb_info(ctxt, me, np)
  if (np == -1) then
    info = psb_err_context_error_
    call psb_errpush(info,name)
    goto 9999
  endif

  ! just this case right now
  ia = 1
  ja = 1

  ix = 1
  if (present(jx)) then
    ijx = jx
  else
    ijx = 1
  endif

  iy = 1
  if (present(jy)) then
    ijy = jy
  else
    ijy = 1
  endif

  if (present(k)) then
    lik = min(k,size(x,2)-ijx+1)
    lik = min(lik,size(y,2)-ijy+1)
  else
    lik = min(size(x,2)-ijx+1,size(y,2)-ijy+1)
  endif

  if (present(choice)) then
    choice_ = choice
  else
    choice_ = psb_avg_
  endif

  if (present(scale)) then
    lscale = psb_toupper(scale)
  else
    lscale = 'U'
  endif

  if (present(trans)) then
    itrans = psb_toupper(trans)
    if((itrans == 'N').or.(itrans == 'T').or.  (itrans == 'C')) then
      ! OK
    else
      info = psb_err_iarg_invalid_value_
      call psb_errpush(info,name)
      goto 9999
    end if
  else
    itrans = 'N'
  endif

  m    = desc_a%get_global_rows()
  nrow = desc_a%get_local_rows()
  ncol = desc_a%get_local_cols()
  lldx = size(x,1)
  lldy = size(y,1)

  if((lldx < ncol).or.(lldy < ncol)) then
    info=psb_err_lld_case_not_implemented_
    call psb_errpush(info,name)
    goto 9999
  end if

  ! check for presence/size of a work area
  iwork => null()
  liwork= 2*ncol
  if (present(work)) then
    if (size(work) >= liwork) then
      aliw =.false.
    else
      aliw=.true.
    endif
  else
    aliw=.true.
  end if

  if (aliw) then
    allocate(iwork(liwork),stat=info)
    if(info /= psb_success_) then
      info=psb_err_from_subroutine_
      ch_err='psb_realloc'
      call psb_errpush(info,name,a_err=ch_err)
      goto 9999
    end if
  else
    iwork => work
  endif

  iwork(1)=0.d0

  if(present(diag)) then
    lld = size(diag)
    id => diag
  else
    lld=1
    allocate(id(1))
    id=1.d0
  end if

  ! checking for matrix correctness
  call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja)
  ! checking for vectors correctness
  if (info == psb_success_) &
       & call psb_chkvect(m,lik,lldx,ix,ijx,desc_a,info,iix,jjx,check_halo=.true.)
  if (info == psb_success_) &
       & call psb_chkvect(m,lik,lldy,iy,ijy,desc_a,info,iiy,jjy,check_halo=.true.)
  if(info /= psb_success_) then
    info=psb_err_from_subroutine_
    ch_err='psb_chkvect/mat'
    call psb_errpush(info,name,a_err=ch_err)
    goto 9999
  end if

  if(ja /= ix) then
    ! this case is not yet implemented
    info = psb_err_ja_nix_ia_niy_unsupported_
  end if

  if((iix /= 1).or.(iiy /= 1)) then
    ! this case is not yet implemented
    info = psb_err_ix_n1_iy_n1_unsupported_
  end if

  if(info /= psb_success_) then
    call psb_errpush(info,name)
    goto 9999
  end if

  ik = lik ! This should not be a problem.
  ! We expect ik to be small, well within IPK
  ! Perform local triangular system solve
  xp => x(iix:lldx,jjx:jjx+ik-1)
  yp => y(iiy:lldy,jjy:jjy+ik-1)
  call psb_cssm(alpha,a,xp,beta,yp,info,scale=scale,d=diag,trans=trans)

  if(info /= psb_success_) then
    info = psb_err_from_subroutine_
    ch_err='cssm'
    call psb_errpush(info,name,a_err=ch_err)
    goto 9999
  end if

  ! update overlap elements
  if (choice_ > 0) then
    call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),ik,&
         & done,yp,desc_a,iwork,info,data=psb_comm_ovr_)

    if (info == psb_success_) call psi_ovrl_upd(yp,desc_a,choice_,info)
    if (info /= psb_success_) then
      call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner updates')
      goto 9999
    end if
  end if

  if(aliw) deallocate(iwork)
  if(.not.present(diag)) deallocate(id)

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(ctxt,err_act)

  return
end subroutine psb_dspsm
!
! Subroutine: psb_dspsv
!  Performs one of the distributed matrix-vector operations
!
!     Y := alpha * Pr * A-1 * Pc * X + beta * Y,   or
!
!     Y := alpha * D * Pr * A-1 * Pc * X + beta * Y, or
!
!     Y := alpha * Pr * A-1 * Pc * D * X + beta * Y,  or
!
!     Y := alpha * Pr * A-T * Pc * X + beta * Y,   or
!
!     Y := alpha * D * Pr * A-T * Pc * X + beta * Y, or
!
!     Y := alpha * Pr * A-T * Pc * D * X + beta * Y,  or
!
!  X is a distributed
!  vector and T is a M-by-M distributed triangular matrix.
!
!
! Arguments:
!    alpha   -  real.               The scalar alpha.
!    a       -  type(psb_dspmat_type). The sparse matrix containing A.
!    x(:)    -  real                The input vector containing the entries of ( X ).
!    beta    -  real                The scalar beta.
!    y(:)    -  real                The input vector containing the entries of ( Y ).
!    desc_a  -  type(psb_desc_type).   The communication descriptor.
!    info    -  integer.               Return code
!    trans   -  character(optional).   Whether A or A'. If not present 'N' is assumed.
!    scale   -  character(optional).   Specify some type of operation with
!                                      the diagonal matrix D.
!    choice  -  integer(optional).     The kind of update to perform on overlap elements.
!    d(:)    -  real, optional      Matrix for diagonal scaling.
!    work(:) -  real, optional      Working area.
!
subroutine  psb_dspsv(alpha,a,x,beta,y,desc_a,info,&
     & trans, scale, choice, diag, work)
  use psb_base_mod, psb_protect_name => psb_dspsv
  use psi_mod
  implicit none

  real(psb_dpk_), intent(in)              :: alpha, beta
  real(psb_dpk_), intent(in), target      :: x(:)
  real(psb_dpk_), intent(inout), target   :: y(:)
  type(psb_dspmat_type), intent(in)         :: a
  type(psb_desc_type), intent(in)           :: desc_a
  integer(psb_ipk_), intent(out)                      :: info
  real(psb_dpk_), intent(in), optional, target    :: diag(:)
  real(psb_dpk_), optional, target, intent(inout) :: work(:)
  character, intent(in), optional           :: trans, scale
  integer(psb_ipk_), intent(in), optional             :: choice

  ! locals
  type(psb_ctxt_type) :: ctxt
  integer(psb_ipk_) :: np, me, &
       & err_act, iix, jjx, iia, jja, lldx,lldy, choice_,&
       & ik, i, lld, nrow, ncol, liwork, llwork, iiy, jjy, idx, ndm
  integer(psb_lpk_) :: ix, ijx, iy, ijy, m, n, ia, ja, lik, jx, jy

  character                :: lscale
  integer(psb_ipk_), parameter       :: nb=4
  real(psb_dpk_),pointer :: iwork(:), xp(:), yp(:), id(:)
  character                :: itrans
  character(len=20)        :: name, ch_err
  logical                  :: aliw

  name='psb_dspsv'
  info=psb_success_
  call psb_erractionsave(err_act)
  if (psb_errstatus_fatal()) then
    info = psb_err_internal_error_ ;    goto 9999
  end if

  ctxt=desc_a%get_context()

  call psb_info(ctxt, me, np)
  if (np == -1) then
    info = psb_err_context_error_
    call psb_errpush(info,name)
    goto 9999
  endif

  ! just this case right now
  ia = 1
  ja = 1
  ix = 1
  iy = 1
  lik = 1
  ik = 1
  jx = 1
  jy = 1

  if (present(choice)) then
    choice_ = choice
  else
    choice_ = psb_avg_
  endif

  if (present(scale)) then
    lscale = psb_toupper(scale)
  else
    lscale = 'U'
  endif

  if (present(trans)) then
    itrans = psb_toupper(trans)
    if((itrans == 'N').or.(itrans == 'T').or.(itrans == 'C')) then
      ! Ok
    else
      info = psb_err_iarg_invalid_value_
      call psb_errpush(info,name)
      goto 9999
    end if
  else
    itrans = 'N'
  endif

  m    = desc_a%get_global_rows()
  nrow = desc_a%get_local_rows()
  ncol = desc_a%get_local_cols()
  lldx = size(x)
  lldy = size(y)

  if((lldx < ncol).or.(lldy < ncol)) then
    info=psb_err_lld_case_not_implemented_
    call psb_errpush(info,name)
    goto 9999
  end if

  iwork => null()
  ! check for presence/size of a work area
  liwork= 2*ncol

  if (present(work)) then
    if (size(work) >= liwork) then
      aliw =.false.
    else
      aliw=.true.
    endif
  else
    aliw=.true.
  end if

  if (aliw) then
    allocate(iwork(liwork),stat=info)
    if(info /= psb_success_) then
      info=psb_err_from_subroutine_
      ch_err='psb_realloc'
      call psb_errpush(info,name,a_err=ch_err)
      goto 9999
    end if
  else
    iwork => work
  endif

  iwork(1)=0.d0

  if(present(diag)) then
    lld = size(diag)
    id => diag
  else
    lld=1
    allocate(id(1))
    id=1.d0
  end if

  ! checking for matrix correctness
  call psb_chkmat(m,m,ia,ja,desc_a,info,iia,jja)
  ! checking for vectors correctness
  if (info == psb_success_) &
       & call psb_chkvect(m,lik,lldx,ix,jx,desc_a,info,iix,jjx,check_halo=.true.)
  if (info == psb_success_) &
       & call psb_chkvect(m,lik,lldy,iy,jy,desc_a,info,iiy,jjy,check_halo=.true.)
  if(info /= psb_success_) then
    info=psb_err_from_subroutine_
    ch_err='psb_chkvect/mat'
    call psb_errpush(info,name,a_err=ch_err)
    goto 9999
  end if

  if(ja /= ix) then
    ! this case is not yet implemented
    info = psb_err_ja_nix_ia_niy_unsupported_
  end if

  if((iix /= 1).or.(iiy /= 1)) then
    ! this case is not yet implemented
    info = psb_err_ix_n1_iy_n1_unsupported_
  end if

  if(info /= psb_success_) then
    call psb_errpush(info,name)
    goto 9999
  end if

  ! Perform local triangular system solve
  xp => x(iix:lldx)
  yp => y(iiy:lldy)
  call psb_cssm(alpha,a,xp,beta,yp,info,scale=scale,d=diag,trans=trans)

  if(info /= psb_success_) then
    info = psb_err_from_subroutine_
    ch_err='dcssm'
    call psb_errpush(info,name,a_err=ch_err)
    goto 9999
  end if

  ! update overlap elements
  if(choice_ > 0) then
    call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),&
         & done,yp,desc_a,iwork,info,data=psb_comm_ovr_)

    if (info == psb_success_) call psi_ovrl_upd(yp,desc_a,choice_,info)
    if (info /= psb_success_) then
      call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner updates')
      goto 9999
    end if
  end if

  if (aliw) deallocate(iwork)
  if(.not.present(diag)) deallocate(id)

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(ctxt,err_act)

  return
end subroutine psb_dspsv


text:
! @@name:	nesting_restrict.4f
! @@type:	F-fixed
! @@compilable:	no
! @@linkable:	no
! @@expect:	failure
      SUBROUTINE WRONG4(N)
      INTEGER N

        INTEGER I
!$OMP   PARALLEL DEFAULT(SHARED)
!$OMP     DO
          DO I = 1, N
            CALL WORK(I, 1)
! incorrect nesting of barrier region in a loop region
!$OMP       BARRIER ! { error "PGF90-S-0155-Illegal context for barrier" }
            CALL WORK(I, 2)
          END DO
!$OMP   END PARALLEL
      END SUBROUTINE WRONG4

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang