Fortran 95 source code to calculate a simple regression curve.(power)

regression_curve_pow.f95

!****************************************************
! 単回帰曲線（べき乗回帰）計算
! : y = a * x**b
! : 連立方程式は ガウスの消去法で解く

!   date          name            version
!   2019.04.10    mk-mode.com     1.00 新規作成
!
! Copyright(C) 2019 mk-mode.com All Rights Reserved.
!****************************************************
!
module const
  ! SP: 単精度(4), DP: 倍精度(8)
  integer,     parameter :: SP = kind(1.0)
  integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
end module const

module comp
  use const
  implicit none
  private
  public :: calc_reg_curve_pow

contains
  ! 単回帰曲線（べき乗回帰）計算
  !
  ! :param(in)  real(8) x(:): 説明変数配列
  ! :param(in)  real(8) y(:): 目的変数配列
  ! :param(out) real(8)    a: 係数 a
  ! :param(out) real(8)    b: 係数 b
  subroutine calc_reg_curve_pow(x, y, a, b)
    implicit none
    real(DP), intent(in)  :: x(:), y(:)
    real(DP), intent(out) :: a, b
    integer(SP) :: size_x, size_y, i
    real(DP)    :: sum_lx, sum_lx2, sum_ly, sum_lxly
    real(DP)    :: mtx(2, 3)
    real(DP), allocatable :: lx(:), ly(:)

    size_x = size(x)
    size_y = size(y)
    if (size_x == 0 .or. size_y == 0) then
      print *, "[ERROR] array size == 0"
      stop
    end if
    if (size_x /= size_y) then
      print *, "[ERROR] size(X) != size(Y)"
      stop
    end if

    allocate(lx(size_x))
    allocate(ly(size_y))
    lx = log(x)
    ly = log(y)
    sum_lx   = sum(lx)
    sum_lx2  = sum(lx * lx)
    sum_ly   = sum(ly)
    sum_lxly = sum(lx * ly)
    deallocate(lx)
    deallocate(ly)
    mtx(1, :) = (/real(size_x, DP),  sum_lx,   sum_ly/)
    mtx(2, :) = (/          sum_lx, sum_lx2, sum_lxly/)
    call solve_ge(2, mtx)
    a = exp(mtx(1, 3))
    b = mtx(2, 3)
  end subroutine calc_reg_curve_pow

  ! 連立方程式を解く（ガウスの消去法）
  !
  ! :param(in)    integer(4)     n: 元数
  ! :param(inout) real(8) a(n,n+1): 係数配列
  subroutine solve_ge(n, a)
    implicit none
    integer(SP), intent(in)    :: n
    real(DP),    intent(inout) :: a(n, n + 1)
    integer(SP) :: i, j
    real(DP)    :: d

    ! 前進消去
    do j = 1, n - 1
      do i = j + 1, n
        d = a(i, j) / a(j, j)
        a(i, j+1:n+1) = a(i, j+1:n+1) - a(j, j+1:n+1) * d
      end do
    end do

    ! 後退代入
    do i = n, 1, -1
      d = a(i, n + 1)
      do j = i + 1, n
        d = d - a(i, j) * a(j, n + 1)
      end do
      a(i, n + 1) = d / a(i, i)
    end do
  end subroutine solve_ge
end module comp

program regression_curve_pow
  use const
  use comp
  implicit none
  character(9), parameter :: F_INP = "input.txt"
  integer(SP),  parameter :: UID   = 10
  real(DP)      :: a, b
  integer(SP)   :: n, i
  character(20) :: f
  real(DP), allocatable :: x(:), y(:)

  ! IN ファイル OPEN
  open(UID, file = F_INP, status = "old")

  ! データ数読み込み
  read(UID, *) n

  ! 配列用メモリ確保
  allocate(x(n))
  allocate(y(n))

  ! データ読み込み
  do i = 1, n
    read(UID, *) x(i), y(i)
  end do
  write(f, '("(A, ", I0, "F8.2, A)")') n
  print f, "説明変数 X = (", x, ")"
  print f, "目的変数 Y = (", y, ")"
  print '(A)', "---"

  ! IN ファイル CLOSE
  close (UID)

  ! 回帰曲線計算
  call calc_reg_curve_pow(x, y, a, b)
  print '(A, F12.8)', "a = ", a
  print '(A, F12.8)', "b = ", b

  ! 配列用メモリ解放
  deallocate(x)
  deallocate(y)

  stop
end program regression_curve_pow