text:
submodule (stdlib_stats) stdlib_stats_corr

  use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
  use stdlib_error, only: error_stop
  use stdlib_linalg, only: diag
  use stdlib_optval, only: optval
  implicit none

contains

    module function corr_1_rsp_rsp(x, dim, mask) result(res)
      real(sp), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      real(sp) :: res

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._sp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_1_rsp_rsp
    module function corr_1_rdp_rdp(x, dim, mask) result(res)
      real(dp), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      real(dp) :: res

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._dp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_1_rdp_rdp
    module function corr_1_rqp_rqp(x, dim, mask) result(res)
      real(qp), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      real(qp) :: res

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._qp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_1_rqp_rqp
    module function corr_1_csp_csp(x, dim, mask) result(res)
      complex(sp), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      real(sp) :: res

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._sp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_1_csp_csp
    module function corr_1_cdp_cdp(x, dim, mask) result(res)
      complex(dp), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      real(dp) :: res

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._dp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_1_cdp_cdp
    module function corr_1_cqp_cqp(x, dim, mask) result(res)
      complex(qp), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      real(qp) :: res

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._qp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_1_cqp_cqp

    module function corr_1_iint8_dp(x, dim, mask) result(res)
      integer(int8), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      real(dp) :: res

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._dp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_1_iint8_dp
    module function corr_1_iint16_dp(x, dim, mask) result(res)
      integer(int16), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      real(dp) :: res

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._dp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_1_iint16_dp
    module function corr_1_iint32_dp(x, dim, mask) result(res)
      integer(int32), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      real(dp) :: res

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._dp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_1_iint32_dp
    module function corr_1_iint64_dp(x, dim, mask) result(res)
      integer(int64), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      real(dp) :: res

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._dp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_1_iint64_dp

    module function corr_mask_1_rsp_rsp(x, dim, mask) result(res)
      real(sp), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:)
      real(sp) :: res

      if (count(mask) < 2) then
        res = ieee_value(1._sp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_mask_1_rsp_rsp
    module function corr_mask_1_rdp_rdp(x, dim, mask) result(res)
      real(dp), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:)
      real(dp) :: res

      if (count(mask) < 2) then
        res = ieee_value(1._dp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_mask_1_rdp_rdp
    module function corr_mask_1_rqp_rqp(x, dim, mask) result(res)
      real(qp), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:)
      real(qp) :: res

      if (count(mask) < 2) then
        res = ieee_value(1._qp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_mask_1_rqp_rqp
    module function corr_mask_1_csp_csp(x, dim, mask) result(res)
      complex(sp), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:)
      real(sp) :: res

      if (count(mask) < 2) then
        res = ieee_value(1._sp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_mask_1_csp_csp
    module function corr_mask_1_cdp_cdp(x, dim, mask) result(res)
      complex(dp), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:)
      real(dp) :: res

      if (count(mask) < 2) then
        res = ieee_value(1._dp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_mask_1_cdp_cdp
    module function corr_mask_1_cqp_cqp(x, dim, mask) result(res)
      complex(qp), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:)
      real(qp) :: res

      if (count(mask) < 2) then
        res = ieee_value(1._qp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_mask_1_cqp_cqp

    module function corr_mask_1_iint8_dp(x, dim, mask) result(res)
      integer(int8), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:)
      real(dp) :: res

      if (count(mask) < 2) then
        res = ieee_value(1._dp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_mask_1_iint8_dp
    module function corr_mask_1_iint16_dp(x, dim, mask) result(res)
      integer(int16), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:)
      real(dp) :: res

      if (count(mask) < 2) then
        res = ieee_value(1._dp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_mask_1_iint16_dp
    module function corr_mask_1_iint32_dp(x, dim, mask) result(res)
      integer(int32), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:)
      real(dp) :: res

      if (count(mask) < 2) then
        res = ieee_value(1._dp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_mask_1_iint32_dp
    module function corr_mask_1_iint64_dp(x, dim, mask) result(res)
      integer(int64), intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:)
      real(dp) :: res

      if (count(mask) < 2) then
        res = ieee_value(1._dp, ieee_quiet_nan)
        return
      end if

      res = 1

    end function corr_mask_1_iint64_dp

    module function corr_2_rsp_rsp(x, dim, mask) result(res)
      real(sp), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      real(sp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                          , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      real(sp) :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
      real(sp) :: center(size(x, 1),size(x, 2))

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._sp, ieee_quiet_nan)
        return
      end if

      mean_ = mean(x, dim)
      select case(dim)
        case(1)
          do i = 1, size(x, 1)
            center(i, :) = x(i, :) - mean_
          end do
            res = matmul( transpose(center), center)
        case(2)
          do i = 1, size(x, 2)
            center(:, i) = x(:, i) - mean_
          end do
            res = matmul( center, transpose(center))
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

      mean_ = 1 / sqrt(diag(res))
      do i = 1, size(res, 1)
        do j = 1, size(res, 2)
          res(j, i) = res(j, i) * mean_(i) * mean_(j)
        end do
      end do

    end function corr_2_rsp_rsp
    module function corr_2_rdp_rdp(x, dim, mask) result(res)
      real(dp), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                          , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      real(dp) :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
      real(dp) :: center(size(x, 1),size(x, 2))

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._dp, ieee_quiet_nan)
        return
      end if

      mean_ = mean(x, dim)
      select case(dim)
        case(1)
          do i = 1, size(x, 1)
            center(i, :) = x(i, :) - mean_
          end do
            res = matmul( transpose(center), center)
        case(2)
          do i = 1, size(x, 2)
            center(:, i) = x(:, i) - mean_
          end do
            res = matmul( center, transpose(center))
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

      mean_ = 1 / sqrt(diag(res))
      do i = 1, size(res, 1)
        do j = 1, size(res, 2)
          res(j, i) = res(j, i) * mean_(i) * mean_(j)
        end do
      end do

    end function corr_2_rdp_rdp
    module function corr_2_rqp_rqp(x, dim, mask) result(res)
      real(qp), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      real(qp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                          , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      real(qp) :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
      real(qp) :: center(size(x, 1),size(x, 2))

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._qp, ieee_quiet_nan)
        return
      end if

      mean_ = mean(x, dim)
      select case(dim)
        case(1)
          do i = 1, size(x, 1)
            center(i, :) = x(i, :) - mean_
          end do
            res = matmul( transpose(center), center)
        case(2)
          do i = 1, size(x, 2)
            center(:, i) = x(:, i) - mean_
          end do
            res = matmul( center, transpose(center))
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

      mean_ = 1 / sqrt(diag(res))
      do i = 1, size(res, 1)
        do j = 1, size(res, 2)
          res(j, i) = res(j, i) * mean_(i) * mean_(j)
        end do
      end do

    end function corr_2_rqp_rqp
    module function corr_2_csp_csp(x, dim, mask) result(res)
      complex(sp), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      complex(sp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                          , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      complex(sp) :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
      complex(sp) :: center(size(x, 1),size(x, 2))

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._sp, ieee_quiet_nan)
        return
      end if

      mean_ = mean(x, dim)
      select case(dim)
        case(1)
          do i = 1, size(x, 1)
            center(i, :) = x(i, :) - mean_
          end do
            res = matmul( transpose(conjg(center)), center)
        case(2)
          do i = 1, size(x, 2)
            center(:, i) = x(:, i) - mean_
          end do
            res = matmul( center, transpose(conjg(center)))
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

      mean_ = 1 / sqrt(diag(res))
      do i = 1, size(res, 1)
        do j = 1, size(res, 2)
          res(j, i) = res(j, i) * mean_(i) * mean_(j)
        end do
      end do

    end function corr_2_csp_csp
    module function corr_2_cdp_cdp(x, dim, mask) result(res)
      complex(dp), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      complex(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                          , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      complex(dp) :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
      complex(dp) :: center(size(x, 1),size(x, 2))

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._dp, ieee_quiet_nan)
        return
      end if

      mean_ = mean(x, dim)
      select case(dim)
        case(1)
          do i = 1, size(x, 1)
            center(i, :) = x(i, :) - mean_
          end do
            res = matmul( transpose(conjg(center)), center)
        case(2)
          do i = 1, size(x, 2)
            center(:, i) = x(:, i) - mean_
          end do
            res = matmul( center, transpose(conjg(center)))
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

      mean_ = 1 / sqrt(diag(res))
      do i = 1, size(res, 1)
        do j = 1, size(res, 2)
          res(j, i) = res(j, i) * mean_(i) * mean_(j)
        end do
      end do

    end function corr_2_cdp_cdp
    module function corr_2_cqp_cqp(x, dim, mask) result(res)
      complex(qp), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      complex(qp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                          , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      complex(qp) :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
      complex(qp) :: center(size(x, 1),size(x, 2))

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._qp, ieee_quiet_nan)
        return
      end if

      mean_ = mean(x, dim)
      select case(dim)
        case(1)
          do i = 1, size(x, 1)
            center(i, :) = x(i, :) - mean_
          end do
            res = matmul( transpose(conjg(center)), center)
        case(2)
          do i = 1, size(x, 2)
            center(:, i) = x(:, i) - mean_
          end do
            res = matmul( center, transpose(conjg(center)))
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

      mean_ = 1 / sqrt(diag(res))
      do i = 1, size(res, 1)
        do j = 1, size(res, 2)
          res(j, i) = res(j, i) * mean_(i) * mean_(j)
        end do
      end do

    end function corr_2_cqp_cqp

    module function corr_2_iint8_dp(x, dim, mask) result(res)
      integer(int8), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                      , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      real(dp) :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
      real(dp) :: center(size(x, 1),size(x, 2))

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._dp, ieee_quiet_nan)
        return
      end if

      mean_ = mean(x, dim)
      select case(dim)
        case(1)
          do i = 1, size(x, 1)
            center(i, :) = real(x(i, :), dp) - mean_
          end do
          res = matmul( transpose(center), center)
        case(2)
          do i = 1, size(x, 2)
            center(:, i) = real(x(:, i), dp) - mean_
          end do
          res = matmul( center, transpose(center))
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

      mean_ = 1 / sqrt(diag(res))
      do i = 1, size(res, 1)
        do j = 1, size(res, 2)
          res(j, i) = res(j, i) * mean_(i) * mean_(j)
        end do
      end do

    end function corr_2_iint8_dp
    module function corr_2_iint16_dp(x, dim, mask) result(res)
      integer(int16), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                      , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      real(dp) :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
      real(dp) :: center(size(x, 1),size(x, 2))

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._dp, ieee_quiet_nan)
        return
      end if

      mean_ = mean(x, dim)
      select case(dim)
        case(1)
          do i = 1, size(x, 1)
            center(i, :) = real(x(i, :), dp) - mean_
          end do
          res = matmul( transpose(center), center)
        case(2)
          do i = 1, size(x, 2)
            center(:, i) = real(x(:, i), dp) - mean_
          end do
          res = matmul( center, transpose(center))
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

      mean_ = 1 / sqrt(diag(res))
      do i = 1, size(res, 1)
        do j = 1, size(res, 2)
          res(j, i) = res(j, i) * mean_(i) * mean_(j)
        end do
      end do

    end function corr_2_iint16_dp
    module function corr_2_iint32_dp(x, dim, mask) result(res)
      integer(int32), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                      , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      real(dp) :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
      real(dp) :: center(size(x, 1),size(x, 2))

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._dp, ieee_quiet_nan)
        return
      end if

      mean_ = mean(x, dim)
      select case(dim)
        case(1)
          do i = 1, size(x, 1)
            center(i, :) = real(x(i, :), dp) - mean_
          end do
          res = matmul( transpose(center), center)
        case(2)
          do i = 1, size(x, 2)
            center(:, i) = real(x(:, i), dp) - mean_
          end do
          res = matmul( center, transpose(center))
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

      mean_ = 1 / sqrt(diag(res))
      do i = 1, size(res, 1)
        do j = 1, size(res, 2)
          res(j, i) = res(j, i) * mean_(i) * mean_(j)
        end do
      end do

    end function corr_2_iint32_dp
    module function corr_2_iint64_dp(x, dim, mask) result(res)
      integer(int64), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                      , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      real(dp) :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
      real(dp) :: center(size(x, 1),size(x, 2))

      if (.not.optval(mask, .true.) .or. size(x) < 2) then
        res = ieee_value(1._dp, ieee_quiet_nan)
        return
      end if

      mean_ = mean(x, dim)
      select case(dim)
        case(1)
          do i = 1, size(x, 1)
            center(i, :) = real(x(i, :), dp) - mean_
          end do
          res = matmul( transpose(center), center)
        case(2)
          do i = 1, size(x, 2)
            center(:, i) = real(x(:, i), dp) - mean_
          end do
          res = matmul( center, transpose(center))
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

      mean_ = 1 / sqrt(diag(res))
      do i = 1, size(res, 1)
        do j = 1, size(res, 2)
          res(j, i) = res(j, i) * mean_(i) * mean_(j)
        end do
      end do

    end function corr_2_iint64_dp

    module function corr_mask_2_rsp_rsp(x, dim, mask) result(res)
      real(sp), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:,:)
      real(sp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                          , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      real(sp) :: centeri_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      real(sp) :: centerj_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      logical :: mask_(merge(size(x, 2), size(x, 1), mask = 1<dim))

      select case(dim)
        case(1)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(:, i) .and. mask(:, j))
             centeri_ = merge( x(:, i) - mean(x(:, i), mask = mask_),&
                0._sp,&
                mask_)
             centerj_ = merge( x(:, j) - mean(x(:, j), mask = mask_),&
                0._sp,&
                mask_)

              res(j, i) = dot_product( centerj_, centeri_)&
               /sqrt(dot_product( centeri_, centeri_)*&
                     dot_product( centerj_, centerj_))

            end do
          end do
        case(2)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(i, :) .and. mask(j, :))
             centeri_ = merge( x(i, :) - mean(x(i, :), mask = mask_),&
                0._sp,&
                mask_)
             centerj_ = merge( x(j, :) - mean(x(j, :), mask = mask_),&
                0._sp,&
                mask_)

              res(j, i) = dot_product( centeri_, centerj_)&
               /sqrt(dot_product( centeri_, centeri_)*&
                     dot_product( centerj_, centerj_))
            end do
          end do
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

    end function corr_mask_2_rsp_rsp
    module function corr_mask_2_rdp_rdp(x, dim, mask) result(res)
      real(dp), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:,:)
      real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                          , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      real(dp) :: centeri_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      real(dp) :: centerj_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      logical :: mask_(merge(size(x, 2), size(x, 1), mask = 1<dim))

      select case(dim)
        case(1)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(:, i) .and. mask(:, j))
             centeri_ = merge( x(:, i) - mean(x(:, i), mask = mask_),&
                0._dp,&
                mask_)
             centerj_ = merge( x(:, j) - mean(x(:, j), mask = mask_),&
                0._dp,&
                mask_)

              res(j, i) = dot_product( centerj_, centeri_)&
               /sqrt(dot_product( centeri_, centeri_)*&
                     dot_product( centerj_, centerj_))

            end do
          end do
        case(2)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(i, :) .and. mask(j, :))
             centeri_ = merge( x(i, :) - mean(x(i, :), mask = mask_),&
                0._dp,&
                mask_)
             centerj_ = merge( x(j, :) - mean(x(j, :), mask = mask_),&
                0._dp,&
                mask_)

              res(j, i) = dot_product( centeri_, centerj_)&
               /sqrt(dot_product( centeri_, centeri_)*&
                     dot_product( centerj_, centerj_))
            end do
          end do
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

    end function corr_mask_2_rdp_rdp
    module function corr_mask_2_rqp_rqp(x, dim, mask) result(res)
      real(qp), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:,:)
      real(qp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                          , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      real(qp) :: centeri_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      real(qp) :: centerj_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      logical :: mask_(merge(size(x, 2), size(x, 1), mask = 1<dim))

      select case(dim)
        case(1)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(:, i) .and. mask(:, j))
             centeri_ = merge( x(:, i) - mean(x(:, i), mask = mask_),&
                0._qp,&
                mask_)
             centerj_ = merge( x(:, j) - mean(x(:, j), mask = mask_),&
                0._qp,&
                mask_)

              res(j, i) = dot_product( centerj_, centeri_)&
               /sqrt(dot_product( centeri_, centeri_)*&
                     dot_product( centerj_, centerj_))

            end do
          end do
        case(2)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(i, :) .and. mask(j, :))
             centeri_ = merge( x(i, :) - mean(x(i, :), mask = mask_),&
                0._qp,&
                mask_)
             centerj_ = merge( x(j, :) - mean(x(j, :), mask = mask_),&
                0._qp,&
                mask_)

              res(j, i) = dot_product( centeri_, centerj_)&
               /sqrt(dot_product( centeri_, centeri_)*&
                     dot_product( centerj_, centerj_))
            end do
          end do
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

    end function corr_mask_2_rqp_rqp
    module function corr_mask_2_csp_csp(x, dim, mask) result(res)
      complex(sp), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:,:)
      complex(sp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                          , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      complex(sp) :: centeri_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      complex(sp) :: centerj_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      logical :: mask_(merge(size(x, 2), size(x, 1), mask = 1<dim))

      select case(dim)
        case(1)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(:, i) .and. mask(:, j))
             centeri_ = merge( x(:, i) - mean(x(:, i), mask = mask_),&
                cmplx(0,0,kind=sp),&
                mask_)
             centerj_ = merge( x(:, j) - mean(x(:, j), mask = mask_),&
                cmplx(0,0,kind=sp),&
                mask_)

              res(j, i) = dot_product( centerj_, centeri_)&
               /sqrt(dot_product( centeri_, centeri_)*&
                     dot_product( centerj_, centerj_))

            end do
          end do
        case(2)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(i, :) .and. mask(j, :))
             centeri_ = merge( x(i, :) - mean(x(i, :), mask = mask_),&
                cmplx(0,0,kind=sp),&
                mask_)
             centerj_ = merge( x(j, :) - mean(x(j, :), mask = mask_),&
                cmplx(0,0,kind=sp),&
                mask_)

              res(j, i) = dot_product( centeri_, centerj_)&
               /sqrt(dot_product( centeri_, centeri_)*&
                     dot_product( centerj_, centerj_))
            end do
          end do
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

    end function corr_mask_2_csp_csp
    module function corr_mask_2_cdp_cdp(x, dim, mask) result(res)
      complex(dp), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:,:)
      complex(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                          , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      complex(dp) :: centeri_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      complex(dp) :: centerj_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      logical :: mask_(merge(size(x, 2), size(x, 1), mask = 1<dim))

      select case(dim)
        case(1)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(:, i) .and. mask(:, j))
             centeri_ = merge( x(:, i) - mean(x(:, i), mask = mask_),&
                cmplx(0,0,kind=dp),&
                mask_)
             centerj_ = merge( x(:, j) - mean(x(:, j), mask = mask_),&
                cmplx(0,0,kind=dp),&
                mask_)

              res(j, i) = dot_product( centerj_, centeri_)&
               /sqrt(dot_product( centeri_, centeri_)*&
                     dot_product( centerj_, centerj_))

            end do
          end do
        case(2)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(i, :) .and. mask(j, :))
             centeri_ = merge( x(i, :) - mean(x(i, :), mask = mask_),&
                cmplx(0,0,kind=dp),&
                mask_)
             centerj_ = merge( x(j, :) - mean(x(j, :), mask = mask_),&
                cmplx(0,0,kind=dp),&
                mask_)

              res(j, i) = dot_product( centeri_, centerj_)&
               /sqrt(dot_product( centeri_, centeri_)*&
                     dot_product( centerj_, centerj_))
            end do
          end do
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

    end function corr_mask_2_cdp_cdp
    module function corr_mask_2_cqp_cqp(x, dim, mask) result(res)
      complex(qp), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:,:)
      complex(qp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                          , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      complex(qp) :: centeri_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      complex(qp) :: centerj_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      logical :: mask_(merge(size(x, 2), size(x, 1), mask = 1<dim))

      select case(dim)
        case(1)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(:, i) .and. mask(:, j))
             centeri_ = merge( x(:, i) - mean(x(:, i), mask = mask_),&
                cmplx(0,0,kind=qp),&
                mask_)
             centerj_ = merge( x(:, j) - mean(x(:, j), mask = mask_),&
                cmplx(0,0,kind=qp),&
                mask_)

              res(j, i) = dot_product( centerj_, centeri_)&
               /sqrt(dot_product( centeri_, centeri_)*&
                     dot_product( centerj_, centerj_))

            end do
          end do
        case(2)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(i, :) .and. mask(j, :))
             centeri_ = merge( x(i, :) - mean(x(i, :), mask = mask_),&
                cmplx(0,0,kind=qp),&
                mask_)
             centerj_ = merge( x(j, :) - mean(x(j, :), mask = mask_),&
                cmplx(0,0,kind=qp),&
                mask_)

              res(j, i) = dot_product( centeri_, centerj_)&
               /sqrt(dot_product( centeri_, centeri_)*&
                     dot_product( centerj_, centerj_))
            end do
          end do
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

    end function corr_mask_2_cqp_cqp

    module function corr_mask_2_iint8_dp(x, dim, mask) result(res)
      integer(int8), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:,:)
      real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                          , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      real(dp) :: centeri_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      real(dp) :: centerj_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      logical :: mask_(merge(size(x, 2), size(x, 1), mask = 1<dim))

      select case(dim)
        case(1)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(:, i) .and. mask(:, j))
             centeri_ = merge( x(:, i) - mean(x(:, i), mask = mask_),&
                0._dp, mask_)
             centerj_ = merge( x(:, j) - mean(x(:, j), mask = mask_),&
                0._dp, mask_)

             res(j, i) = dot_product( centerj_, centeri_)&
                 /sqrt(dot_product( centeri_, centeri_)*&
                       dot_product( centerj_, centerj_))

            end do
          end do
        case(2)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(i, :) .and. mask(j, :))
             centeri_ = merge( x(i, :) - mean(x(i, :), mask = mask_),&
                0._dp, mask_)
             centerj_ = merge( x(j, :) - mean(x(j, :), mask = mask_),&
                0._dp, mask_)

             res(j, i) = dot_product( centeri_, centerj_)&
                 /sqrt(dot_product( centeri_, centeri_)*&
                       dot_product( centerj_, centerj_))
            end do
          end do
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

    end function corr_mask_2_iint8_dp
    module function corr_mask_2_iint16_dp(x, dim, mask) result(res)
      integer(int16), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:,:)
      real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                          , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      real(dp) :: centeri_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      real(dp) :: centerj_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      logical :: mask_(merge(size(x, 2), size(x, 1), mask = 1<dim))

      select case(dim)
        case(1)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(:, i) .and. mask(:, j))
             centeri_ = merge( x(:, i) - mean(x(:, i), mask = mask_),&
                0._dp, mask_)
             centerj_ = merge( x(:, j) - mean(x(:, j), mask = mask_),&
                0._dp, mask_)

             res(j, i) = dot_product( centerj_, centeri_)&
                 /sqrt(dot_product( centeri_, centeri_)*&
                       dot_product( centerj_, centerj_))

            end do
          end do
        case(2)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(i, :) .and. mask(j, :))
             centeri_ = merge( x(i, :) - mean(x(i, :), mask = mask_),&
                0._dp, mask_)
             centerj_ = merge( x(j, :) - mean(x(j, :), mask = mask_),&
                0._dp, mask_)

             res(j, i) = dot_product( centeri_, centerj_)&
                 /sqrt(dot_product( centeri_, centeri_)*&
                       dot_product( centerj_, centerj_))
            end do
          end do
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

    end function corr_mask_2_iint16_dp
    module function corr_mask_2_iint32_dp(x, dim, mask) result(res)
      integer(int32), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:,:)
      real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                          , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      real(dp) :: centeri_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      real(dp) :: centerj_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      logical :: mask_(merge(size(x, 2), size(x, 1), mask = 1<dim))

      select case(dim)
        case(1)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(:, i) .and. mask(:, j))
             centeri_ = merge( x(:, i) - mean(x(:, i), mask = mask_),&
                0._dp, mask_)
             centerj_ = merge( x(:, j) - mean(x(:, j), mask = mask_),&
                0._dp, mask_)

             res(j, i) = dot_product( centerj_, centeri_)&
                 /sqrt(dot_product( centeri_, centeri_)*&
                       dot_product( centerj_, centerj_))

            end do
          end do
        case(2)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(i, :) .and. mask(j, :))
             centeri_ = merge( x(i, :) - mean(x(i, :), mask = mask_),&
                0._dp, mask_)
             centerj_ = merge( x(j, :) - mean(x(j, :), mask = mask_),&
                0._dp, mask_)

             res(j, i) = dot_product( centeri_, centerj_)&
                 /sqrt(dot_product( centeri_, centeri_)*&
                       dot_product( centerj_, centerj_))
            end do
          end do
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

    end function corr_mask_2_iint32_dp
    module function corr_mask_2_iint64_dp(x, dim, mask) result(res)
      integer(int64), intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:,:)
      real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                          , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i, j
      real(dp) :: centeri_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      real(dp) :: centerj_(merge(size(x, 2), size(x, 1), mask = 1<dim))
      logical :: mask_(merge(size(x, 2), size(x, 1), mask = 1<dim))

      select case(dim)
        case(1)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(:, i) .and. mask(:, j))
             centeri_ = merge( x(:, i) - mean(x(:, i), mask = mask_),&
                0._dp, mask_)
             centerj_ = merge( x(:, j) - mean(x(:, j), mask = mask_),&
                0._dp, mask_)

             res(j, i) = dot_product( centerj_, centeri_)&
                 /sqrt(dot_product( centeri_, centeri_)*&
                       dot_product( centerj_, centerj_))

            end do
          end do
        case(2)
          do i = 1, size(res, 2)
            do j = 1, size(res, 1)
             mask_ = merge(.true., .false., mask(i, :) .and. mask(j, :))
             centeri_ = merge( x(i, :) - mean(x(i, :), mask = mask_),&
                0._dp, mask_)
             centerj_ = merge( x(j, :) - mean(x(j, :), mask = mask_),&
                0._dp, mask_)

             res(j, i) = dot_product( centeri_, centerj_)&
                 /sqrt(dot_product( centeri_, centeri_)*&
                       dot_product( centerj_, centerj_))
            end do
          end do
        case default
          call error_stop("ERROR (corr): wrong dimension")
      end select

    end function corr_mask_2_iint64_dp

end submodule

max_stars_repo_name: LKedward/stdlib-fpm

max_issues_repo_name: LKedward/stdlib-fpm

max_forks_repo_name: LKedward/stdlib-fpm


text:
program H2O_vacf
   implicit none
   integer,parameter:: nstep=40000
   character :: b
   real(8) :: velocity(nstep,216*9),cor_dip(150000),Tmax,dt
   real(8) :: sumu,T=298.15d0,vacf
   integer :: i,j,k ,N,tau
   integer,parameter:: tau_nstep=4000
   real(8)::elec(3*216),vtmp(3),dotdip(tau_nstep,3),uacf(0:tau_nstep)
   open(14,file='mdvel')
   read(14,*) b
   do j=1,nstep
      read(14,*) velocity(j,:)
   enddo
   print *,'finish read velocity'

   do i=0,215
      elec(3*i+1)=-2.d0
      elec(3*i+2)=1.d0
      elec(3*i+3)=1.d0
   enddo
   ! calculate dot (dip) auto-correlation functions
   ! uacf(0:ncorr), vel(0:nsteps, ndof)
   call dotdipacf(uacf, velocity, elec, 216*3, 9, nstep, tau_nstep)

   open(13,file='vacf_test.txt')
   do j=0,tau_nstep
      write(13,*) uacf(j)
   enddo
   close(13)

end program

subroutine dotdipacf(uacf, vel, elec, natoms, ndof, nsteps, ncorr)
   integer :: ncorr
   integer :: natoms
   integer :: ndof
   real(8) :: vel(nsteps, ndof)
   real(8) :: elec(natoms)
   real(8) :: uacf(0:ncorr)
   real(8) :: acf(0:ncorr, 3)

   integer :: i, j, k
   real(8) :: dotdip(nsteps,3)
   real(8) :: vtmp(3)

   ! loop over all steps
   do i = 1, nsteps
      vtmp = 0d0
      k = 0
      ! loop over all atoms
      do j = 1, natoms
         vtmp(1) = vtmp(1) + elec(j) * vel(i, k+1)
         vtmp(2) = vtmp(2) + elec(j) * vel(i, k+2)
         vtmp(3) = vtmp(3) + elec(j) * vel(i, k+3)
         k = k + 3
      end do
      dotdip(i,1:3) = vtmp(1:3)
   end do

   do i = 1,3
      call autocorr(acf(0:ncorr,i), dotdip(1:nsteps,i), nsteps, ncorr)
   end do
   do i = 0, ncorr
      uacf(i) = acf(i,1) + acf(i,2) + acf(i,3)
   end do

   return
end subroutine dotdipacf

subroutine autocorr(acf, A, iTtotal, iTcorr)
   implicit none
   integer, intent(in)  ::  iTtotal       ! total number of A
   integer, intent(in)  ::  iTcorr        ! correlation length
   real(8), intent(in)  ::  A(1:iTtotal)  ! array of values of A
   real(8), intent(out) ::  acf(0:iTcorr) ! normalized auto-correlation function

   real(8) ::  nACF(0:iTcorr)             ! number of ACF, for normalization of ACF
   integer :: iT, iT0, iT0max, iTau       ! time index of t0+tau, t0, t0max and tau

   ! initialize all the auto-correlation functions to be calculated
   acf(0:iTcorr)  = 0d0
   nACF(0:iTcorr) = 0d0

   ! for each t0 from 1 to total, calculate auto-correlation function
   do iT0 = 1, iTtotal

      ! the max index of T0 must not exceed the last index
      iT0max = min(iTtotal, iT0+iTcorr)

      ! for all possible t0+tau, calculate a(t0)*a(t0+tau) and add them
      do iT = iT0, iT0max

         ! interval between t and t0
         iTau  = iT - iT0

         ! calculate auto-correlation function for current interval
         acf(iTau) = acf(iTau) + A(iT0) * A(iT)

         ! count the number of auto-correlation function
         nACF(iTau) = nACF(iTau) + 1d0
      end do

   end do

   ! normalize auto-correlation function
   do iT = 0, iTcorr
      acf(iT) = acf(iT) / nACF(iT)
   end do

   return
end subroutine autocorr

max_stars_repo_name: jianliugroup/LSC-IVR-Amber

max_issues_repo_name: jianliugroup/LSC-IVR-Amber

max_forks_repo_name: jianliugroup/LSC-IVR-Amber


text:
program test_prime
!! test correctness of primality detection
use, intrinsic :: iso_fortran_env, only: int64, dp=>real64
use numerical, only: isprime
implicit none

integer, allocatable :: M(:)
integer :: N, i, j
integer(int64) :: tic, toc, rate
character(10) :: argv

call system_clock(count_rate=rate)

N = 100
call get_command_argument(1,argv, status=i)
if(i==0) read(argv,*) N

allocate(M(N))
j=1
i=2
print *,'finding first ',N, 'primes'

call system_clock(tic)
do while(j<=N)
  if (isprime(i)) then
    M(j) = i
    j = j + 1
  endif
  i = i + 1
enddo
call system_clock(toc)

print '(A,ES10.2,A)',' in ',(toc-tic) / real(rate,dp),' seconds.'
print '(10000I12)', M

!! TODO show logrithmic growth rate (increasing prime sparsity)

end program

max_stars_repo_name: scivision/preschool-coding

max_issues_repo_name: scivision/preschool-coding

max_forks_repo_name: scivision/preschool-coding


text:
C$Procedure DAFRWA ( DAF, record/word to address )

      SUBROUTINE DAFRWA ( RECNO, WORDNO, ADDR )

C$ Abstract
C
C     Convert a record/word pair to its equivalent address within
C     a DAF.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     DAF
C
C$ Keywords
C
C     CONVERSION
C     FILES
C
C$ Declarations

      INTEGER               RECNO
      INTEGER               WORDNO
      INTEGER               ADDR

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     RECNO,
C     WORDNO     I   Record, word numbers of a location within DAF.
C     ADDR       O   Corresponding address.
C
C$ Detailed_Input
C
C     RECNO,
C     WORDNO      are the record and word numbers of an arbitrary
C                 location within a DAF.
C
C$ Detailed_Output
C
C     ADDR        is the corresponding address within the DAF.
C
C$ Parameters
C
C      None.
C
C$ Files
C
C     None.
C
C$ Exceptions
C
C     1) If either RECNO or WORDNO is zero or negative, the error
C        SPICE(DAFNOSUCHADDR) is signalled.
C
C$ Particulars
C
C     To the user, the data in a DAF appear to be a contiguous
C     collection of double precision numbers, each of which has an
C     address. To the DAF software, however, the data appear to be
C     a collection of records, each containing 128 double precision
C     words. The routines DAFARW and DAFRWA translate between these
C     two representations.
C
C$ Examples
C
C     Routines DAFRDA and DAFWDA illustrate the use of DAFARW and
C     DAFRWA.
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     NAIF Document 167.0, "Double Precision Array Files (DAF)
C     Specification and User's Guide"
C
C$ Author_and_Institution
C
C     I.M. Underwood  (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.0.2, 10-MAR-1992 (WLT)
C
C        Comment section for permuted index source lines was added
C        following the header.
C
C-    SPICELIB Version 1.0.1, 22-MAR-1990 (HAN)
C
C        Literature references added to the header.
C
C-    SPICELIB Version 1.0.0, 31-JAN-1990 (IMU)
C
C-&

C$ Index_Entries
C
C     record/word to daf address
C
C-&

C
C     SPICELIB functions
C
      LOGICAL               RETURN

C
C     Standard SPICE error handling.
C
      IF ( RETURN () ) THEN
         RETURN

      ELSE IF ( RECNO .LE. 0  .OR.  WORDNO .LE. 0 ) THEN
         CALL CHKIN  ( 'DAFRWA'                           )
         CALL SETMSG ( 'No address for record #, word #.' )
         CALL ERRINT ( '#', RECNO                         )
         CALL ERRINT ( '#', WORDNO                        )
         CALL SIGERR ( 'SPICE(DAFNOSUCHADDR)'             )
         CALL CHKOUT ( 'DAFRWA'                           )
         RETURN
      END IF

C
C     If the record and word numbers are legal, the computation is
C     straightforward.
C
      ADDR = WORDNO + ( RECNO - 1 ) * 128

      RETURN

C$Procedure DAFARW ( DAF, address to record/word )

      ENTRY DAFARW ( ADDR, RECNO, WORDNO )

C$ Abstract
C
C     Convert an address within a DAF to its equivalent
C     record/word representation.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     DAF
C
C$ Keywords
C
C     CONVERSION
C     FILES
C
C$ Declarations
C
C     INTEGER               ADDR
C     INTEGER               RECNO
C     INTEGER               WORDNO
C
C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     ADDR       I   Address within DAF.
C     RECNO,
C     WORDNO     O   Corresponding record, word numbers.
C
C$ Detailed_Input
C
C     ADDR        is an arbitrary address within a DAF.
C
C$ Detailed_Output
C
C     RECNO,
C     WORDNO      are the corresponding record and word numbers
C                 within the DAF.
C
C$ Parameters
C
C      None.
C
C$ Files
C
C     None.
C
C$ Exceptions
C
C     1) If ADDR is zero or negative, the error SPICE(DAFNOSUCHADDR)
C        is signalled.
C
C$ Particulars
C
C     To the user, the data in a DAF appear to be a contiguous
C     collection of double precision numbers, each of which has an
C     address. To the DAF software, however, the data appear to be
C     a collection of records, each containing 128 double precision
C     words. The routines DAFARW and DAFRWA translate between these
C     two representations.
C
C$ Examples
C
C     Routines DAFRDA and DAFWDA illustrate the use of DAFARW and
C     DAFRWA.
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     NAIF Document 167.0, "Double Precision Array Files (DAF)
C     Specification and User's Guide"
C
C$ Author_and_Institution
C
C     I.M. Underwood  (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.0.2, 10-MAR-1992 (WLT)
C
C        Comment section for permuted index source lines was added
C        following the header.
C
C-    SPICELIB Version 1.0.1, 22-MAR-1990 (HAN)
C
C        Literature references added to the header.
C
C-    SPICELIB Version 1.0.0, 31-JAN-1990 (IMU)
C
C-&

C$ Index_Entries
C
C     daf address to record/word
C
C-&

C
C     Standard SPICE error handling.
C
      IF ( RETURN () ) THEN
         RETURN

      ELSE IF ( ADDR .LE. 0 ) THEN
         CALL CHKIN  ( 'DAFARW'                         )
         CALL SETMSG ( 'No record, word for address #.' )
         CALL ERRINT ( '#', ADDR                        )
         CALL SIGERR ( 'SPICE(DAFNOSUCHADDR)'           )
         CALL CHKOUT ( 'DAFARW'                         )
         RETURN
      END IF

C
C     If the address is legal, the computation is straightforward.
C
      RECNO  =        ( ADDR  - 1 ) / 128 + 1
      WORDNO = ADDR - ( RECNO - 1 ) * 128

      RETURN
      END

ext: f

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
*DECK LA05AS
      SUBROUTINE LA05AS (A, IND, NZ, IA, N, IP, IW, W, G, U)
C***BEGIN PROLOGUE  LA05AS
C***SUBSIDIARY
C***PURPOSE  Subsidiary to SPLP
C***LIBRARY   SLATEC
C***TYPE      SINGLE PRECISION (LA05AS-S, LA05AD-D)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C     THIS SUBPROGRAM IS A SLIGHT MODIFICATION OF A SUBPROGRAM
C     FROM THE C. 1979 AERE HARWELL LIBRARY.  THE NAME OF THE
C     CORRESPONDING HARWELL CODE CAN BE OBTAINED BY DELETING
C     THE FINAL LETTER =S= IN THE NAMES USED HERE.
C     REVISIONS MADE BY R J HANSON, SNLA, AUGUST, 1979.
C     REVISED SEP. 13, 1979.
C
C     ROYALTIES HAVE BEEN PAID TO AERE-UK FOR USE OF THEIR CODES
C     IN THE PACKAGE GIVEN HERE.  ANY PRIMARY USAGE OF THE HARWELL
C     SUBROUTINES REQUIRES A ROYALTY AGREEMENT AND PAYMENT BETWEEN
C     THE USER AND AERE-UK.  ANY USAGE OF THE SANDIA WRITTEN CODES
C     SPLP( ) (WHICH USES THE HARWELL SUBROUTINES) IS PERMITTED.
C
C IP(I,1),IP(I,2) POINT TO THE START OF ROW/COL I.
C IW(I,1),IW(I,2) HOLD THE NUMBER OF NON-ZEROS IN ROW/COL I.
C DURING THE MAIN BODY OF THIS SUBROUTINE THE VECTORS IW(.,3),IW(.,5),
C     IW(.,7) ARE USED TO HOLD DOUBLY LINKED LISTS OF ROWS THAT HAVE
C     NOT BEEN PIVOTAL AND HAVE EQUAL NUMBERS OF NON-ZEROS.
C IW(.,4),IW(.,6),IW(.,8) HOLD SIMILAR LISTS FOR THE COLUMNS.
C IW(I,3),IW(I,4) HOLD FIRST ROW/COLUMN TO HAVE I NON-ZEROS
C     OR ZERO IF THERE ARE NONE.
C IW(I,5), IW(I,6) HOLD ROW/COL NUMBER OF ROW/COL PRIOR TO ROW/COL I
C     IN ITS LIST, OR ZERO IF NONE.
C IW(I,7), IW(I,8) HOLD ROW/COL NUMBER OF ROW/COL AFTER ROW/COL I
C     IN ITS LIST, OR ZERO IF NONE.
C FOR ROWS/COLS THAT HAVE BEEN PIVOTAL IW(I,5),IW(I,6) HOLD NEGATION OF
C     POSITION OF ROW/COL I IN THE PIVOTAL ORDERING.
C
C***SEE ALSO  SPLP
C***ROUTINES CALLED  LA05ES, MC20AS, R1MACH, XERMSG, XSETUN
C***COMMON BLOCKS    LA05DS
C***REVISION HISTORY  (YYMMDD)
C   811215  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890605  Corrected references to XERRWV.  (WRB)
C   890831  Modified array declarations.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900402  Added TYPE section.  (WRB)
C   900510  Convert XERRWV calls to XERMSG calls.  (RWC)
C***END PROLOGUE  LA05AS
      INTEGER IP(N,2)
      INTEGER IND(IA,2), IW(N,8)
      REAL A(*), AMAX, AU, AM, G, U, SMALL, W(*)
      LOGICAL FIRST
      CHARACTER*8 XERN0, XERN1, XERN2
C
      COMMON /LA05DS/ SMALL, LP, LENL, LENU, NCP, LROW, LCOL
C EPS IS THE RELATIVE ACCURACY OF FLOATING-POINT COMPUTATION
      SAVE EPS, FIRST
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  LA05AS
      IF (FIRST) THEN
         EPS = 2.0E0 * R1MACH(4)
      ENDIF
      FIRST = .FALSE.
C
C     SET THE OUTPUT UNIT NUMBER FOR THE ERROR PROCESSOR.
C     THE USAGE OF THIS ERROR PROCESSOR IS DOCUMENTED IN THE
C     SANDIA LABS. TECH. REPT. SAND78-1189, BY R E JONES.
      CALL XSETUN(LP)
      IF (U.GT.1.0E0) U = 1.0E0
      IF (U.LT.EPS) U = EPS
      IF (N.LT.1) GO TO 670
      G = 0.
      DO 50 I=1,N
         W(I) = 0.
         DO 40 J=1,5
            IW(I,J) = 0
   40    CONTINUE
   50 CONTINUE
C
C FLUSH OUT SMALL ENTRIES, COUNT ELEMENTS IN ROWS AND COLUMNS
      L = 1
      LENU = NZ
      DO 80 IDUMMY=1,NZ
         IF (L.GT.LENU) GO TO 90
         DO 60 K=L,LENU
            IF (ABS(A(K)).LE.SMALL) GO TO 70
            I = IND(K,1)
            J = IND(K,2)
            G = MAX(ABS(A(K)),G)
            IF (I.LT.1 .OR. I.GT.N) GO TO 680
            IF (J.LT.1 .OR. J.GT.N) GO TO 680
            IW(I,1) = IW(I,1) + 1
            IW(J,2) = IW(J,2) + 1
   60    CONTINUE
         GO TO 90
   70    L = K
         A(L) = A(LENU)
         IND(L,1) = IND(LENU,1)
         IND(L,2) = IND(LENU,2)
         LENU = LENU - 1
   80 CONTINUE
C
   90 LENL = 0
      LROW = LENU
      LCOL = LROW
C MCP IS THE MAXIMUM NUMBER OF COMPRESSES PERMITTED BEFORE AN
C     ERROR RETURN RESULTS.
      MCP = MAX(N/10,20)
      NCP = 0
C CHECK FOR NULL ROW OR COLUMN AND INITIALIZE IP(I,2) TO POINT
C     JUST BEYOND WHERE THE LAST COMPONENT OF COLUMN I OF A WILL
C     BE STORED.
      K = 1
      DO 110 IR=1,N
         K = K + IW(IR,2)
         IP(IR,2) = K
         DO 100 L=1,2
            IF (IW(IR,L).LE.0) GO TO 700
  100    CONTINUE
  110 CONTINUE
C REORDER BY ROWS
C CHECK FOR DOUBLE ENTRIES WHILE USING THE NEWLY CONSTRUCTED
C     ROW FILE TO CONSTRUCT THE COLUMN FILE. NOTE THAT BY PUTTING
C    THE ENTRIES IN BACKWARDS AND DECREASING IP(J,2) EACH TIME IT
C     IS USED WE AUTOMATICALLY LEAVE IT POINTING TO THE FIRST ELEMENT.
      CALL MC20AS(N, LENU, A, IND(1,2), IP, IND(1,1), 0)
      KL = LENU
      DO 130 II=1,N
         IR = N + 1 - II
         KP = IP(IR,1)
         DO 120 K=KP,KL
            J = IND(K,2)
            IF (IW(J,5).EQ.IR) GO TO 660
            IW(J,5) = IR
            KR = IP(J,2) - 1
            IP(J,2) = KR
            IND(KR,1) = IR
  120    CONTINUE
         KL = KP - 1
  130 CONTINUE
C
C SET UP LINKED LISTS OF ROWS AND COLS WITH EQUAL NUMBERS OF NON-ZEROS.
      DO 150 L=1,2
         DO 140 I=1,N
            NZ = IW(I,L)
            IN = IW(NZ,L+2)
            IW(NZ,L+2) = I
            IW(I,L+6) = IN
            IW(I,L+4) = 0
            IF (IN.NE.0) IW(IN,L+4) = I
  140    CONTINUE
  150 CONTINUE
C
C
C START OF MAIN ELIMINATION LOOP.
      DO 590 IPV=1,N
C FIND PIVOT. JCOST IS MARKOWITZ COST OF CHEAPEST PIVOT FOUND SO FAR,
C     WHICH IS IN ROW IPP AND COLUMN JP.
         JCOST = N*N
C LOOP ON LENGTH OF COLUMN TO BE SEARCHED
         DO 240 NZ=1,N
            IF (JCOST.LE.(NZ-1)**2) GO TO 250
            J = IW(NZ,4)
C SEARCH COLUMNS WITH NZ NON-ZEROS.
            DO 190 IDUMMY=1,N
               IF (J.LE.0) GO TO 200
               KP = IP(J,2)
               KL = KP + IW(J,2) - 1
               DO 180 K=KP,KL
                  I = IND(K,1)
                  KCOST = (NZ-1)*(IW(I,1)-1)
                  IF (KCOST.GE.JCOST) GO TO 180
                  IF (NZ.EQ.1) GO TO 170
C FIND LARGEST ELEMENT IN ROW OF POTENTIAL PIVOT.
                  AMAX = 0.
                  K1 = IP(I,1)
                  K2 = IW(I,1) + K1 - 1
                  DO 160 KK=K1,K2
                     AMAX = MAX(AMAX,ABS(A(KK)))
                     IF (IND(KK,2).EQ.J) KJ = KK
  160             CONTINUE
C PERFORM STABILITY TEST.
                  IF (ABS(A(KJ)).LT.AMAX*U) GO TO 180
  170             JCOST = KCOST
                  IPP = I
                  JP = J
                  IF (JCOST.LE.(NZ-1)**2) GO TO 250
  180          CONTINUE
               J = IW(J,8)
  190       CONTINUE
C SEARCH ROWS WITH NZ NON-ZEROS.
  200       I = IW(NZ,3)
            DO 230 IDUMMY=1,N
               IF (I.LE.0) GO TO 240
               AMAX = 0.
               KP = IP(I,1)
               KL = KP + IW(I,1) - 1
C FIND LARGEST ELEMENT IN THE ROW
               DO 210 K=KP,KL
                  AMAX = MAX(ABS(A(K)),AMAX)
  210          CONTINUE
               AU = AMAX*U
               DO 220 K=KP,KL
C PERFORM STABILITY TEST.
                  IF (ABS(A(K)).LT.AU) GO TO 220
                  J = IND(K,2)
                  KCOST = (NZ-1)*(IW(J,2)-1)
                  IF (KCOST.GE.JCOST) GO TO 220
                  JCOST = KCOST
                  IPP = I
                  JP = J
                  IF (JCOST.LE.(NZ-1)**2) GO TO 250
  220          CONTINUE
               I = IW(I,7)
  230       CONTINUE
  240    CONTINUE
C
C PIVOT FOUND.
C REMOVE ROWS AND COLUMNS INVOLVED IN ELIMINATION FROM ORDERING VECTORS.
  250    KP = IP(JP,2)
         KL = IW(JP,2) + KP - 1
         DO 290 L=1,2
            DO 280 K=KP,KL
               I = IND(K,L)
               IL = IW(I,L+4)
               IN = IW(I,L+6)
               IF (IL.EQ.0) GO TO 260
               IW(IL,L+6) = IN
               GO TO 270
  260          NZ = IW(I,L)
               IW(NZ,L+2) = IN
  270          IF (IN.GT.0) IW(IN,L+4) = IL
  280       CONTINUE
            KP = IP(IPP,1)
            KL = KP + IW(IPP,1) - 1
  290    CONTINUE
C STORE PIVOT
         IW(IPP,5) = -IPV
         IW(JP,6) = -IPV
C ELIMINATE PIVOTAL ROW FROM COLUMN FILE AND FIND PIVOT IN ROW FILE.
         DO 320 K=KP,KL
            J = IND(K,2)
            KPC = IP(J,2)
            IW(J,2) = IW(J,2) - 1
            KLC = KPC + IW(J,2)
            DO 300 KC=KPC,KLC
               IF (IPP.EQ.IND(KC,1)) GO TO 310
  300       CONTINUE
  310       IND(KC,1) = IND(KLC,1)
            IND(KLC,1) = 0
            IF (J.EQ.JP) KR = K
  320    CONTINUE
C BRING PIVOT TO FRONT OF PIVOTAL ROW.
         AU = A(KR)
         A(KR) = A(KP)
         A(KP) = AU
         IND(KR,2) = IND(KP,2)
         IND(KP,2) = JP
C
C PERFORM ELIMINATION ITSELF, LOOPING ON NON-ZEROS IN PIVOT COLUMN.
         NZC = IW(JP,2)
         IF (NZC.EQ.0) GO TO 550
         DO 540 NC=1,NZC
            KC = IP(JP,2) + NC - 1
            IR = IND(KC,1)
C SEARCH NON-PIVOT ROW FOR ELEMENT TO BE ELIMINATED.
            KR = IP(IR,1)
            KRL = KR + IW(IR,1) - 1
            DO 330 KNP=KR,KRL
               IF (JP.EQ.IND(KNP,2)) GO TO 340
  330       CONTINUE
C BRING ELEMENT TO BE ELIMINATED TO FRONT OF ITS ROW.
  340       AM = A(KNP)
            A(KNP) = A(KR)
            A(KR) = AM
            IND(KNP,2) = IND(KR,2)
            IND(KR,2) = JP
            AM = -A(KR)/A(KP)
C COMPRESS ROW FILE UNLESS IT IS CERTAIN THAT THERE IS ROOM FOR NEW ROW.
            IF (LROW+IW(IR,1)+IW(IPP,1)+LENL.LE.IA) GO TO 350
            IF (NCP.GE.MCP .OR. LENU+IW(IR,1)+IW(IPP,1)+LENL.GT.IA) GO
     *       TO 710
            CALL LA05ES(A, IND(1,2), IP, N, IW, IA, .TRUE.)
            KP = IP(IPP,1)
            KR = IP(IR,1)
  350       KRL = KR + IW(IR,1) - 1
            KQ = KP + 1
            KPL = KP + IW(IPP,1) - 1
C PLACE PIVOT ROW (EXCLUDING PIVOT ITSELF) IN W.
            IF (KQ.GT.KPL) GO TO 370
            DO 360 K=KQ,KPL
               J = IND(K,2)
               W(J) = A(K)
  360       CONTINUE
  370       IP(IR,1) = LROW + 1
C
C TRANSFER MODIFIED ELEMENTS.
            IND(KR,2) = 0
            KR = KR + 1
            IF (KR.GT.KRL) GO TO 430
            DO 420 KS=KR,KRL
               J = IND(KS,2)
               AU = A(KS) + AM*W(J)
               IND(KS,2) = 0
C IF ELEMENT IS VERY SMALL REMOVE IT FROM U.
               IF (ABS(AU).LE.SMALL) GO TO 380
               G = MAX(G,ABS(AU))
               LROW = LROW + 1
               A(LROW) = AU
               IND(LROW,2) = J
               GO TO 410
  380          LENU = LENU - 1
C REMOVE ELEMENT FROM COL FILE.
               K = IP(J,2)
               KL = K + IW(J,2) - 1
               IW(J,2) = KL - K
               DO 390 KK=K,KL
                  IF (IND(KK,1).EQ.IR) GO TO 400
  390          CONTINUE
  400          IND(KK,1) = IND(KL,1)
               IND(KL,1) = 0
  410          W(J) = 0.
  420       CONTINUE
C
C SCAN PIVOT ROW FOR FILLS.
  430       IF (KQ.GT.KPL) GO TO 520
            DO 510 KS=KQ,KPL
               J = IND(KS,2)
               AU = AM*W(J)
               IF (ABS(AU).LE.SMALL) GO TO 500
               LROW = LROW + 1
               A(LROW) = AU
               IND(LROW,2) = J
               LENU = LENU + 1
C
C CREATE FILL IN COLUMN FILE.
               NZ = IW(J,2)
               K = IP(J,2)
               KL = K + NZ - 1
               IF (NZ .EQ. 0) GO TO 460
C IF POSSIBLE PLACE NEW ELEMENT AT END OF PRESENT ENTRY.
               IF (KL.NE.LCOL) GO TO 440
               IF (LCOL+LENL.GE.IA) GO TO 460
               LCOL = LCOL + 1
               GO TO 450
  440          IF (IND(KL+1,1).NE.0) GO TO 460
  450          IND(KL+1,1) = IR
               GO TO 490
C NEW ENTRY HAS TO BE CREATED.
  460          IF (LCOL+LENL+NZ+1.LT.IA) GO TO 470
C COMPRESS COLUMN FILE IF THERE IS NOT ROOM FOR NEW ENTRY.
               IF (NCP.GE.MCP .OR. LENU+LENL+NZ+1.GE.IA) GO TO 710
               CALL LA05ES(A, IND, IP(1,2), N, IW(1,2), IA, .FALSE.)
               K = IP(J,2)
               KL = K + NZ - 1
C TRANSFER OLD ENTRY INTO NEW.
  470          IP(J,2) = LCOL + 1
               IF (KL .LT. K) GO TO 485
               DO 480 KK=K,KL
                  LCOL = LCOL + 1
                  IND(LCOL,1) = IND(KK,1)
                  IND(KK,1) = 0
  480          CONTINUE
  485          CONTINUE
C ADD NEW ELEMENT.
               LCOL = LCOL + 1
               IND(LCOL,1) = IR
  490          G = MAX(G,ABS(AU))
               IW(J,2) = NZ + 1
  500          W(J) = 0.
  510       CONTINUE
  520       IW(IR,1) = LROW + 1 - IP(IR,1)
C
C STORE MULTIPLIER
            IF (LENL+LCOL+1.LE.IA) GO TO 530
C COMPRESS COL FILE IF NECESSARY.
            IF (NCP.GE.MCP) GO TO 710
            CALL LA05ES(A, IND, IP(1,2), N, IW(1,2), IA, .FALSE.)
  530       K = IA - LENL
            LENL = LENL + 1
            A(K) = AM
            IND(K,1) = IPP
            IND(K,2) = IR
            LENU = LENU - 1
  540    CONTINUE
C
C INSERT ROWS AND COLUMNS INVOLVED IN ELIMINATION IN LINKED LISTS
C     OF EQUAL NUMBERS OF NON-ZEROS.
  550    K1 = IP(JP,2)
         K2 = IW(JP,2) + K1 - 1
         IW(JP,2) = 0
         DO 580 L=1,2
            IF (K2.LT.K1) GO TO 570
            DO 560 K=K1,K2
               IR = IND(K,L)
               IF (L.EQ.1) IND(K,L) = 0
               NZ = IW(IR,L)
               IF (NZ.LE.0) GO TO 720
               IN = IW(NZ,L+2)
               IW(IR,L+6) = IN
               IW(IR,L+4) = 0
               IW(NZ,L+2) = IR
               IF (IN.NE.0) IW(IN,L+4) = IR
  560       CONTINUE
  570       K1 = IP(IPP,1) + 1
            K2 = IW(IPP,1) + K1 - 2
  580    CONTINUE
  590 CONTINUE
C
C RESET COLUMN FILE TO REFER TO U AND STORE ROW/COL NUMBERS IN
C     PIVOTAL ORDER IN IW(.,3),IW(.,4)
      DO 600 I=1,N
         J = -IW(I,5)
         IW(J,3) = I
         J = -IW(I,6)
         IW(J,4) = I
         IW(I,2) = 0
  600 CONTINUE
      DO 620 I=1,N
         KP = IP(I,1)
         KL = IW(I,1) + KP - 1
         DO 610 K=KP,KL
            J = IND(K,2)
            IW(J,2) = IW(J,2) + 1
  610    CONTINUE
  620 CONTINUE
      K = 1
      DO 630 I=1,N
         K = K + IW(I,2)
         IP(I,2) = K
  630 CONTINUE
      LCOL = K - 1
      DO 650 II=1,N
         I = IW(II,3)
         KP = IP(I,1)
         KL = IW(I,1) + KP - 1
         DO 640 K=KP,KL
            J = IND(K,2)
            KN = IP(J,2) - 1
            IP(J,2) = KN
            IND(KN,1) = I
  640    CONTINUE
  650 CONTINUE
      RETURN
C
C     THE FOLLOWING INSTRUCTIONS IMPLEMENT THE FAILURE EXITS.
C
  660 IF (LP.GT.0) THEN
         WRITE (XERN1, '(I8)') IR
         WRITE (XERN2, '(I8)') J
         CALL XERMSG ('SLATEC', 'LA05AS', 'MORE THAN ONE MATRIX ' //
     *      'ENTRY.  HERE ROW = ' // XERN1 // ' AND COL = ' // XERN2,
     *      -4, 1)
      ENDIF
      G = -4.
      RETURN
C
  670 IF (LP.GT.0) CALL XERMSG ('SLATEC', 'LA05AS',
     *   'THE ORDER OF THE SYSTEM, N, IS NOT POSITIVE.', -1, 1)
      G = -1.0E0
      RETURN
C
  680 IF (LP.GT.0) THEN
         WRITE (XERN0, '(I8)') K
         WRITE (XERN1, '(I8)') I
         WRITE (XERN2, '(I8)') J
         CALL XERMSG ('SLATEC', 'LA05AS', 'ELEMENT K = ' // XERN0 //
     *      ' IS OUT OF BOUNDS.$$HERE ROW = ' // XERN1 //
     *      ' AND COL = ' // XERN2, -3, 1)
      ENDIF
      G = -3.
      RETURN
C
  700 IF (LP.GT.0) THEN
         WRITE (XERN1, '(I8)') L
         CALL XERMSG ('SLATEC', 'LA05AS', 'ROW OR COLUMN HAS NO ' //
     *      'ELEMENTS.  HERE INDEX = ' // XERN1, -2, 1)
      ENDIF
      G = -2.
      RETURN
C
  710 IF (LP.GT.0) CALL XERMSG ('SLATEC', 'LA05AS',
     *   'LENGTHS OF ARRAYS A(*) AND IND(*,2) ARE TOO SMALL.', -7, 1)
      G = -7.
      RETURN
C
  720 IPV = IPV + 1
      IW(IPV,1) = IR
      DO 730 I=1,N
         II = -IW(I,L+4)
         IF (II.GT.0) IW(II,1) = I
  730 CONTINUE
C
      IF (LP.GT.0) THEN
         XERN1 = 'ROWS'
         IF (L.EQ.2) XERN1 = 'COLUMNS'
         CALL XERMSG ('SLATEC', 'LA05AS', 'DEPENDANT ' // XERN1, -5, 1)
C
  740    WRITE (XERN1, '(I8)') IW(I,1)
         XERN2 = ' '
         IF (I+1.LE.IPV) WRITE (XERN2, '(I8)') IW(I+1,1)
         CALL XERMSG ('SLATEC', 'LA05AS',
     *      'DEPENDENT VECTOR INDICES ARE ' // XERN1 // ' AND ' //
     *      XERN2, -5, 1)
         I = I + 2
         IF (I.LE.IPV) GO TO 740
      ENDIF
      G = -5.
      RETURN
      END


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      subroutine set_effi( temper, radiat, ext   , depthw, daylen,
     &                     id    )
!>\file
!>       calculate and store efficiency for all species

!     modules

      USE      DATA_3DL   ! data and routine for 3D light approach

!     BLOOM commons

      INCLUDE 'blmdim.inc'
      INCLUDE 'size.inc'
      INCLUDE 'phyt2.inc'
      INCLUDE 'putin1.inc'
      INCLUDE 'arran.inc'

!     arguments

      REAL     TEMPER     ! input , temperature
      REAL     RADIAT     ! input , radiation
      REAL     EXT        ! input , total extinction
      REAL     DEPTHW     ! input , depth of the layer
      REAL     DAYLEN     ! input , daylength in hours
      INTEGER  ID         ! input , weeknumber

!     local decalarations

      REAL*8   ALPHA      ! reflection factor
      REAL*8   TEMP       ! temperature
      REAL*8   CSOL       ! radiation
      REAL*8   DSOL       ! radiation
      REAL*8   DEP        ! depth
      REAL*8   EXTTOT     ! total extinction
      REAL*8   DAY        ! daylength in hours
      REAL*8   DEAT       ! DEAT
      REAL*8   ZOODD      ! ZOODD
      REAL*8   TCORR      ! TCORR
      REAL*8   SURF_TYP   ! scaled, converted and corrected radiation for a type
      INTEGER  IGROUP     ! index number of BLOOM algae group
      INTEGER  ITYPE      ! index number of BLOOM algae type
      REAL*8   PMAX20(MT),SDMIXN(MT)

      DEP    = DEPTHW
      EXTTOT = EXT
      TEMP   = TEMPER
      CSOL   = SOLACO * RADIAT
      DAY    = DAYLEN
      DEAT   = 0D0
      ZOODD  = 0D0
      CALL NATMOR ( DEAT  , ZOODD , TEMP  , 1)
      DO ITYPE = 1,NTYP_3DL
         IF (SDMIX(ITYPE) .LT. 0.0) THEN
            SDMIXN(ITYPE) = 1.0D0 + SDMIX(ITYPE)
!           DMIX(K) = DABS(SDMIX(ITYPE)) * DEP
         ELSE
            SDMIXN(ITYPE) = 0.0D0
         ENDIF
      ENDDO

      CALL MAXPRD ( TEFCUR )
      DO ITYPE = 1,NTYP_3DL
         PMAX20(ITYPE) = PMAX(ITYPE)
      ENDDO
      CALL MAXPRD ( TEMP  )

      DSOL=1428.57*CSOL
      DO IGROUP = 1 , NGRO_3DL
         DO ITYPE = IT2(IGROUP,1),IT2(IGROUP,2)
            TCORR      = PMAX20(ITYPE)/PMAX(ITYPE)
            SURF_TYP   = TCORR * DSOL * DEXP (- EXTTOT * SDMIXN(ITYPE) * DEP)
            SURF_TYP   = SURF_TYP/DAY
            CALL EFFILAY_3DL( SURF_TYP, EXTTOT, DEP   , IGROUP, ITYPE )
         ENDDO
      ENDDO

      RETURN
      END SUBROUTINE SET_EFFI

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
! This Source Code Form is subject to the terms of the Mozilla Public
! License, v. 2.0. If a copy of the MPL was not distributed with this
! file, You can obtain one at http://mozilla.org/MPL/2.0/.
module iaib_mod
use IOBuffer
contains
  subroutine make_iaib(cho_data)
    ! A subroutine to construct the (ia|ib) integrals for the sigmavector routines

    use global_var_mod
    use molecule_var_mod
    use cholesky_structs
    use utilities_mod
    use wp_tov_mod
#ifdef TIGER_USE_OMP
    use omp_lib
#endif

    implicit none

    type(cholesky_data),intent(in)::cho_data

    integer::i,a,b,ia_count,status,offset,ab_count,iaib_count
    real(real8),dimension(:,:),allocatable::ia_cho
    real(real8)::int2

    real(kind=real8), external::ddot ! BLAS1 (dot product)

#ifdef TIGER_USE_OMP
    integer::threadID
    threadID = OMP_get_thread_num()+1
#else
    integer,parameter::threadID = 1
#endif

    allocate(ia_cho(numcho,num_orbitals-num_internal),stat=status)
    call allocatecheck(status,"make_iaib")

    iaib_count = 0
    do i = 1, num_internal

       ia_cho = 0.0D0
       do a = 1, num_external

          if (ignorable_pair(a+num_internal,i) ) cycle

          ia_count = a+num_internal
          ia_count = ia_count*(ia_count-1)/2+i
          ia_count = cho_data%mo_ind_inv(ia_count)

          if ( ia_count .eq. 0 ) cycle
          call for_double_buf_readblock(mo_int_no, ia_count, ia_cho(1:numcho,a), threadID)
       enddo

       offset = num_internal*num_external*(num_external+1)/2 + (i-1)*num_external*(num_external+1)/2
       ab_count = 0

       do a = 1, num_external
          do b = 1,a

             ab_count = ab_count + 1

             if (ignorable_pair(a+num_internal,i) ) cycle
             if (ignorable_pair(b+num_internal,i) ) cycle

             iaib_count = iaib_count + 1
             int2 = ddot(numcho,ia_cho(1,a),1,ia_cho(1,b),1)

             call for_double_buf_writeElement(cd_iaib_no,iaib_count,int2,threadID)

          enddo
       enddo
    enddo

    deallocate(ia_cho,stat=status)
    call deallocatecheck(status,"make_iaib")

  end subroutine make_iaib
end module iaib_mod

max_stars_repo_name: EACcodes/TigerCI

max_issues_repo_name: EACcodes/TigerCI

max_forks_repo_name: EACcodes/TigerCI


text:
C$Procedure   STCG01 ( STAR catalog type 1, get star data )

      SUBROUTINE STCG01 ( INDEX,  RA,     DEC,    RASIG,
     .                    DECSIG, CATNUM, SPTYPE, VMAG )

C$ Abstract
C
C     Get data for a single star from a SPICE type 1 star catalog.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     EK
C
C$ Keywords
C
C     None.
C
C$ Declarations

      INTEGER               INDEX
      DOUBLE PRECISION      RA
      DOUBLE PRECISION      DEC
      DOUBLE PRECISION      RASIG
      DOUBLE PRECISION      DECSIG
      INTEGER               CATNUM
      CHARACTER*(*)         SPTYPE
      DOUBLE PRECISION      VMAG

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     INDEX       I   Star index.
C     RA          O   Right ascension in radians.
C     DEC         O   Declination in radians.
C     RAS         O   Right ascension uncertainty in radians.
C     DECS        O   Declination uncertainty in radians.
C     CATNUM      O   Catalog number.
C     SPTYPE      O   Spectral type.
C     VMAG        O   Visual magnitude.
C
C$ Detailed_Input
C
C     INDEX       is the index of the star in the list of stars
C                 that satisfy the selection criteria specified in
C                 the last call to STCF01.
C
C$ Detailed_Output
C
C     RA          is right ascension of the star at the catalog epoch
C                 in radians relative to the J2000 inertial frame.
C
C     DEC         is declination of the star at the catalog epoch in
C                 radians relative to the J2000 inertial frame.
C
C     RASIG       is the uncertainty in right ascension of the star at
C                 the catalog epoch in radians.
C
C     DECSIG      is the uncertainty in declination of the star at
C                 the catalog epoch in radians.
C
C     CATNUM      is the star number in the catalog.
C
C     SPTYPE      is the star's spectral type. See catalog description
C                 for more information regarding encoding of spectral
C                 type values.
C
C     VMAG        is the visual magnitude of the star.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1) If fetching of any of output values fails, then
C        the error 'SPICE(BADSTARINDEX)' is signalled.
C
C     2) If no star catalog has been loaded, the error is dianosed
C        by a routine called by this one.
C
C     3) If STCF01 was not called first, the EK query
C        error 'SPICE(INVALIDINDEX)' is signalled.
C
C$ Files
C
C     This routine reads the data from SPICE type 1 star catalog file
C     loaded into the program by a call to STCL01.
C
C     SPICE type 1 star catalog files MUST contain a single data table.
C     It can occupy a single segment or it can spread across multiple
C     segments. This table MUST include the following columns:
C
C        column name                data type          units
C     -------------------------------------------------------
C        RA                   DOUBLE PRECISION        DEGREES
C        DEC                  DOUBLE PRECISION        DEGREES
C        RA_SIGMA             DOUBLE PRECISION        DEGREES
C        DEC_SIGMA            DOUBLE PRECISION        DEGREES
C        CATALOG_NUMBER       INTEGER
C        SPECTRAL_TYPE        CHARACTER*(4)
C        VISUAL_MAGNITUDE     DOUBLE PRECISION
C
C     Nulls are not allowed in any of the columns.
C     Other columns can also be present in the table but their data
C     will NOT be accessible through STCF01 and STCG01 --
C     the interface used to access data in the catalog. Note
C     that the names and attributes of these additional columns
C     must be identical for all segments containing this table.
C
C$ Particulars
C
C     This routine is intended to be a part of the user interface to
C     the SPICE type 1 star catalog. It allows the caller to retrieve
C     data for a single star found by STCF01 using the star's
C     index within the search result array. This subroutine MUST
C     NOT be called before a search by STCF01 was done.
C
C     Other routines in the SPICE type 1 star catalog access
C     family are:
C
C        STCL01  load the catalog file and make its data
C                available for search and retrieval.
C
C        STCF01  search through the catalog for all stars within
C                a specified RA-DEC rectangle.
C
C$ Examples
C
C     In the following code fragment, STCG01 is used to retrieve
C     position and characteristics for every star within an RA - DEC
C     rectangle from a particular SPICE type 1 star catalog.
C
C     C
C     C     Load catalog file.
C     C
C           CALL STCL01 ( CATFN, TABNAM, HANDLE )
C     C
C     C     Search through the loaded catalog.
C     C
C           CALL STCF01 ( TABNAM, RAMIN,  RAMAX,
C          .              DECMIN, DECMAX, NSTARS )
C     C
C     C     Retrieve data for every star found.
C     C
C           DO I = 1, NSTARS
C
C              CALL STCG01 ( I, RA, DEC, RASIG, DECSIG,
C          .                 CATNUM, SPTYPE, VMAG )
C
C           END DO
C
C$ Restrictions
C
C     1) The catalog file STCG01 reads data from MUST be loaded
C        by STCL01 and a search through the catalog MUST be done by
C        STCF01 before STCG01 is called.
C
C     2) No other EK queries can be made between the call to STCF01
C        and the call to STCG01.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     B.V. Semenov    (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.0.0, 15-MAY-1996 (BVS)
C
C-&

C$ Index_Entries
C
C     get data for single star from a type 1 star catalog
C
C-&

C$ Revisions
C
C-&

C
C
C     SPICELIB functions
C
      LOGICAL               RETURN
      DOUBLE PRECISION      RPD

C
C     Local variables.
C
      LOGICAL               FOUND
      LOGICAL               NULL

C
C     Standard SPICE error handling.
C
      IF ( RETURN () ) THEN
         RETURN
      ELSE
         CALL CHKIN ( 'STCG01' )
      END IF

C
C     Fetch data from the catalog in the following order
C     as defined QUERY string template in STCF01 routine
C
C           RA, DEC, RASIG, DECSIG, CATNUM, SPTYPE, VMAG
C
C     Check FOUNDs and report error if any of the parameters
C     is not found.
C
C     Since NULLs are not allowed in any of the star catalog
C     columns, no check for NULLs is performed.
C
      CALL EKGD (  1, INDEX, 1, RA, NULL, FOUND )
      IF ( .NOT. FOUND ) THEN
         CALL SETMSG ( 'RA value for star # not found. '    )
         CALL ERRINT ( '#', INDEX                           )
         CALL SIGERR ( 'SPICE(BADSTARINDEX)'                )
         CALL CHKOUT ( 'STCG01' )
         RETURN
      END IF

      CALL EKGD (  2, INDEX, 1, DEC, NULL, FOUND )
      IF ( .NOT. FOUND ) THEN
         CALL SETMSG ( 'DEC value for star # not found. '   )
         CALL ERRINT ( '#', INDEX                           )
         CALL SIGERR ( 'SPICE(BADSTARINDEX)'                )
         CALL CHKOUT ( 'STCG01' )
         RETURN
      END IF

      CALL EKGD (  3, INDEX, 1, RASIG, NULL, FOUND )
      IF ( .NOT. FOUND ) THEN
         CALL SETMSG ( 'RASIG value for star # not found. ' )
         CALL ERRINT ( '#', INDEX                           )
         CALL SIGERR ( 'SPICE(BADSTARINDEX)'                )
         CALL CHKOUT ( 'STCG01' )
         RETURN
      END IF

      CALL EKGD (  4, INDEX, 1, DECSIG, NULL, FOUND )
      IF ( .NOT. FOUND ) THEN
         CALL SETMSG ( 'DECSIG value for star # not found.' )
         CALL ERRINT ( '#', INDEX                           )
         CALL SIGERR ( 'SPICE(BADSTARINDEX)'                )
         CALL CHKOUT ( 'STCG01' )
         RETURN
      END IF

      CALL EKGI (  5, INDEX, 1, CATNUM, NULL, FOUND )
      IF ( .NOT. FOUND ) THEN
         CALL SETMSG ( 'CATNUM value for star # not found.' )
         CALL ERRINT ( '#', INDEX                           )
         CALL SIGERR ( 'SPICE(BADSTARINDEX)'                )
         CALL CHKOUT ( 'STCG01' )
         RETURN
      END IF

      CALL EKGC (  6, INDEX, 1, SPTYPE, NULL, FOUND )
      IF ( .NOT. FOUND ) THEN
         CALL SETMSG ( 'SPTYPE value for star # not found.' )
         CALL ERRINT ( '#', INDEX                           )
         CALL SIGERR ( 'SPICE(BADSTARINDEX)'                )
         CALL CHKOUT ( 'STCG01' )
         RETURN
      END IF

      CALL EKGD (  7, INDEX, 1, VMAG, NULL, FOUND )
      IF ( .NOT. FOUND ) THEN
         CALL SETMSG ( 'VMAG value for star # not found. '  )
         CALL ERRINT ( '#', INDEX                           )
         CALL SIGERR ( 'SPICE(BADSTARINDEX)'                )
         CALL CHKOUT ( 'STCG01' )
         RETURN
      END IF

C
C     Convert angles to radians before return.
C
      RA     = RA     * RPD()
      DEC    = DEC    * RPD()
      RASIG  = RASIG  * RPD()
      DECSIG = DECSIG * RPD()

      CALL CHKOUT ( 'STCG01' )
      RETURN
      END


text:
module quadratures
  !-----------------------------------------------------------------------------
  ! Module: quadratures
  !  This modules contains subroutines for computing and testing numerical
  !  quadratures.
  !
  ! Subroutines:
  !  gausslegendre (public)
  !   Compute the Gauss-Legendre quadrature of degree N using the Golub-Welsh
  !   algorithm.
  !
  !  gausslegendre_shifted (public)
  !   Compute the Gauss-Legendre quadrature of degree N using the Golub-Welsh
  !   algorithm and map to the interval (0, 1).
  !
  !  gausslegendre_truncated (public)
  !   Compute the Gauss-Legendre quadrature of degree 2N using the Golub-Welsh
  !   algorithm and truncate returning only the positive gridpoints.
  !
  !  polytest (public)
  !   Test an integral quadrature against the analytic solution of an integral
  !   of the form integral[a,b] x^{m} dx
  !
  ! Useful Citations:
  !  Golub, Gene H., and John H. Welsch. "Calculation of Gauss
  !  quadrature rules." Mathematics of computation (1969): 221-230.
  !---
  ! Developed by J.P. Dark
  ! Contact: email@jpdark.com
  ! Last modified June 2017
  !-----------------------------------------------------------------------------

  implicit none

  integer, parameter :: sp = kind(1.0)
  integer, parameter :: dp = kind(1.d0)
  real(dp), parameter :: PI = 3.141592653589793_dp

  private
  public :: gausslegendre, gausslegendre_truncated, gausslegendre_shifted, &
       clencurt, polytest

contains

  subroutine gausslegendre(N, mu, wt)
    !!------------------------------------------------------------------------
    !! Name: gausslegendre
    !!  Compute the Gauss-Legendre quadrature of degree N using the
    !!  Golub-Welsh algorithm.
    !!
    !! Synopsis
    !!  subroutine gausslegendre(N, mu, wt)
    !!
    !!    INPUT
    !!     Integer ------------- N
    !!
    !!    OUTPUT
    !!     Double precision ---- mu(N), wt(N)
    !!
    !! Purpose
    !!  Compute the Nth order Gauss-Legendre with quadrature points mu and
    !!  the associated weigths wt on the interval (-1, 1).
    !!
    !!  Citation:
    !!  Golub, Gene H., and John H. Welsch. "Calculation of Gauss
    !!  quadrature rules." Mathematics of computation (1969): 221-230.
    !!
    !!
    !! Arguments N (input) integer The number of quadrature points.
    !!
    !!  mu (output) double precision, array(N)
    !!    Quadrature points.
    !!
    !!  wt (output) double precision, array(N)
    !!    Associated weights.
    !!------------------------------------------------------------------------
    implicit none
    integer(8), intent(in) :: N
    real(dp), intent(out), dimension(N) :: mu, wt
    real(dp) :: z(N, N)
    real(dp) :: e(N-1)
    real(dp) :: work(N-1, N-1)
    real(dp) :: tmp
    integer(8) :: info
    integer(8) :: mm
    if (mod(N, 2) /= 0) then
       write(*, *) &
            '> WARNING: Quadrature contains the origin.'
    endif
    mu = 0.0_dp
    do mm = 1, N-1
       tmp = 1.0_dp / (2.0_dp * real(mm, dp)) ** 2
       e(mm) = 0.5_dp / dsqrt(1.0_dp - tmp)
    enddo
    call dsteqr('I', N, mu, e, z, N, work, info)
    do mm = 1, N
       wt(mm) = 2.0_dp * z(1, mm) ** 2
    enddo
  end subroutine gausslegendre

  subroutine gausslegendre_shifted(N, mu, wt)
    !!------------------------------------------------------------------------
    !! Name: gausslegendre_shifted
    !!  Compute the Gauss-Legendre quadrature of degree N using the
    !!  Golub-Welsh algorithm and map to the interval (0, 1)
    !!
    !! Synopsis
    !!  subroutine gausslegendre(N, mu, wt)
    !!
    !!    INPUT
    !!     Integer ------------- N
    !!
    !!    OUTPUT
    !!     Double precision ---- mu(N), wt(N)
    !!
    !! Purpose
    !!  Compute the Nth order Gauss-Legendre with quadrature points mu and
    !!  the associated weigths wt on the interval (-1, 1) and shift using
    !!  linear transform
    !!
    !!     mu   -->  0.5 * (mu + 1)
    !!
    !!  into the interval (0, 1).
    !!
    !!  Citation:
    !!   Golub, Gene H., and John H. Welsch. "Calculation of Gauss
    !!   quadrature rules." Mathematics of computation (1969): 221-230.
    !!
    !!
    !! Arguments
    !!     N (input) integer
    !!       The number of quadrature points.
    !!
    !!     mu (output) double precision, array(N)
    !!       Quadrature points.
    !!
    !!     wt (output) double precision, array(N)
    !!       Associated weights.
    !!------------------------------------------------------------------------
    implicit none
    integer(8), intent(in) :: N
    real(dp), intent(out) :: mu(N), wt(N)
    real(dp) :: z(N, N)
    real(dp) :: e(N-1)
    real(dp) :: work(N-1, N-1)
    real(dp) :: tmp
    integer(8) :: info
    integer(8) :: mm
    mu = 0.0_dp
    do mm = 1, N-1
       tmp = 1.0_dp / (2.0_dp * real(mm, dp)) ** 2
       e(mm) = 0.5_dp / dsqrt(1.0_dp - tmp)
    enddo
    call dsteqr('I', N, mu, e, z, N, work, info)
    mu = 0.5_dp * (mu + 1.0_dp)
    do mm = 1, N
       wt(mm) = z(1, mm) ** 2
    enddo
  end subroutine gausslegendre_shifted

  subroutine gausslegendre_truncated(N, mu, wt)
    !!--------------------------------------------------------------------------
    !! Name: gausslegendre_truncated
    !!  Compute the Gauss-Legendre quadrature of degree 2N using the
    !!  Golub-Welsh algorithm and truncate returning the positive gridpoints.
    !!
    !! Synopsis
    !!  subroutine gausslegendre(N, mu, wt)
    !!
    !!    INPUT
    !!     Integer ------------- N
    !!
    !!    OUTPUT
    !!     Double precision ---- mu(N), wt(N)
    !!
    !! Purpose
    !!  Compute the (2*N)th order Gauss-Legendre with quadrature points mu and
    !!  the associated weigths wt on the interval (-1, 1) and return the N
    !!  gridpoints and associated weights on the interval (0, 1).
    !!
    !!     Citation:
    !!     Golub, Gene H., and John H. Welsch. "Calculation of Gauss
    !!     quadrature rules." Mathematics of computation (1969): 221-230.
    !!
    !!
    !! Arguments
    !!  N (input) integer
    !!    The number of quadrature points.
    !!
    !!  mu (output) double precision, array(N)
    !!    Quadrature points.
    !!
    !!  wt (output) double precision, array(N)
    !!    Associated weights.
    !!---------------------------------------------------------------------------
    implicit none
    integer(8), intent(in) :: N
    real(dp), dimension(N) :: mu, wt
    real(dp), dimension(2*N) :: x
    real(dp) :: z(2*N, 2*N)
    real(dp) :: e(2*N-1)
    real(dp) :: work(2*N-1, 2*N-1)
    real(dp) :: tmp
    integer(8) :: info
    integer(8) :: mm
    x = 0.0_dp
    do mm = 1, 2*N-1
       tmp = 1.0_dp / (2.0_dp * real(mm, dp)) ** 2
       e(mm) = 0.5_dp / dsqrt(1.0_dp - tmp)
    enddo
    call dsteqr('I', 2*N, x, e, z, 2*N, work, info)
    mu = x(N+1:2*N)
    do mm = 1, N
       wt(mm) = 2.0_dp * z(1, N + mm) ** 2
    enddo
  end subroutine gausslegendre_truncated

  subroutine clencurt(N, mu, wt)
    !!--------------------------------------------------------------------------
    !! Name: gausslegendre_truncated
    !!  Compute the Clenshaw-Curtis quadrature of degree N for even N.
    !!
    !! Synopsis
    !!  subroutine clencurt(N, mu, wt)
    !!
    !!    INPUT
    !!     Integer ------------- N
    !!
    !!    OUTPUT
    !!     Double precision ---- mu(N), wt(N)
    !!
    !! Purpose
    !!  Compute the Clenshaw-Curtis quadrature on the interval [-1, 1].
    !!
    !!   *Method adapted from Matlab function in book:
    !!       Spectral Methods in MATLAB by Lloyd N. Trefethen (p.128).
    !!
    !! Arguments
    !!  N (input) integer
    !!    The number of quadrature points (even).
    !!
    !!  mu (output) double precision, array(N)
    !!    Quadrature points.
    !!
    !!  wt (output) double precision, array(N)
    !!    Associated weights.
    !!---------------------------------------------------------------------------
    implicit none
    integer(8), intent(in) :: N
    real(dp), intent(out), dimension(N) :: mu, wt
    real(dp), dimension(N) :: theta
    integer(8) :: M
    real(dp) :: tmp1
    integer(8) :: ii, jj

    !---------------------------------------------------------------------------
    ! Chebyshev gridpoints mu = cos(theta) where theta[n] = pi * (n-1) / N.
    !
    theta(1) = 0.0_dp
    tmp1 = PI / real(N-1, dp)
    do ii = 1, N-1
       theta(ii+1) = tmp1 * real(ii, dp)
    enddo
    mu = cos(theta)
    !---------------------------------------------------------------------------

    !---------------------------------------------------------------------------
    ! Compute weight
    !
    M = N-1
    if (mod(N, 2) .eq. 0) then
       wt(1) = 1.0_dp / real((N-1) ** 2, dp)
       wt(N) = wt(1)
       wt(2:M) = 1.0_dp
       do jj = 1, N/2-1
          ! wt -= 2 * cos(2 * j * theta) / (4 * j^2 - 1)
          tmp1 = real(2 * jj ** 2, dp) - 0.5_dp
          wt(2:M) = wt(2:M) - cos(2.0_dp * real(jj, dp) * theta(2:M)) &
               / tmp1
       enddo
       ! wt *= 2 / (M)
       wt(2:M) = wt(2:M) / (real(N/2, dp) - 0.5_dp)
    else
       write(*, '(A)') '> Warning in __quadratures_MOD_clencurt'
       write(*, '(A)') &
            '> Warning: Odd number of quadrature points - origin is included.'

       wt(1) = 1.0_dp / real(M ** 2 - 1, dp)
       wt(N) = wt(1)
       wt(2:M) = 1.0_dp
       do jj = 1, M/2 - 1
          ! wt -= 2 * cos(2 * j * theta) / (4 * j^2 - 1)
          tmp1 = real(2 * jj ** 2, dp) - 0.5_dp
          wt(2:M) = wt(2:M) - cos(2.0_dp * real(jj, dp) * theta(2:M)) &
               / tmp1
       enddo
       ! wt -= cos((M) * theta) / ((M)^2 - 1)
       wt(2:M) = wt(2:M) - cos(real(M, dp) * theta(2:M)) &
            / real(M ** 2 - 1, dp)
       ! wt *= 2 / (M)
       wt(2:M) = wt(2:M) / real(M/2, dp)
    endif

  end subroutine clencurt

  subroutine polytest(N, m, a, b, x, wt)
    !!--------------------------------------------------------------------------
    !! Name: polytest
    !!  Test an integral quadrature against the analytic solution of an integral
    !!  of the form integral[a,b] x^{m} dx
    !!
    !!
    !! Synopsis
    !!     subroutine polytest( N, M, a, b, x, wt )
    !!
    !!     Integer ------------- N, M
    !!     Double precision ---- a, b, x( N ), wt( N )
    !!
    !! Purpose
    !!  Compute the integral
    !!
    !!     I = Integral_{a}^{b} x^{M} dx
    !!
    !!  both analytically,
    !!
    !!     I_true = 1/(M+1) (b^{M+1} - a^{M+1}),
    !!
    !!  and with the quadrature,
    !!
    !!     I_quad = sum[n=1,N] x(n)^m * wt(n).
    !!
    !!  Compare the results to verify the correctness of the quadrature.
    !!
    !! Arguments
    !!  N (input) integer
    !!    The number of quadrature points.
    !!
    !!  M (input) integer
    !!    Order of the polynomial in integral.
    !!
    !!  a (output) double precision, array(N)
    !!    Start point of the interval the quadrature is defined on.
    !!
    !!  b (output) double precision, array(N)
    !!    End point of the interval the quadrature is defined on.
    !!
    !!  x (output) double precision, array(N)
    !!    Quadrature points.
    !!
    !!  wt (output) double precision, array(N)
    !!    Associated weights.
    !!---------------------------------------------------------------------------
    implicit none
    ! Input
    integer(8), intent(in) :: N
    integer, intent(in) :: M
    real(dp), intent(in) :: a, b
    real(dp), intent(in), dimension(N) :: x, wt
    ! Workspace
    real(dp) :: approximate, analytic
    if (M .eq. -1) then
       analytic = log(b) - log(a)
    else
       analytic = (b ** (M+1) - a ** (M+1)) / real(M+1, dp)
    endif
    approximate = sum(x ** M * wt)
    write(*, '(A, F4.1, A, F4.1, A, I4)') &
         '> Integrate x^m from b to a where: a = ', a, ', b = ',  b, &
         ', m =  ', M
    write(*, '(A, F8.6,A, F8.6, A, E12.4)') 'Analytic = ', analytic, &
    '; Numeric =  ', approximate, '; Error = ', analytic - approximate
    write(*,*) ""
  end subroutine polytest

end module quadratures

max_stars_repo_name: ugks/pddom

max_issues_repo_name: ugks/pddom

max_forks_repo_name: ugks/pddom