dc_units.f90

Go to the documentation of this file.

!== dc_units.f90 - 単位系処理用モジュール
!
! Authors::   Eizi TOYODA, Yasuhiro MORIKAWA
! Version::   $Id: dc_units.f90,v 1.1 2009-03-20 09:09:52 morikawa Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2000-2005. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
! This file provides dc_units
!

module dc_units !:nodoc:
  !
  !== Overview
  !
  ! このモジュールは単位系処理のためのプログラム群を提供します。
  !

  use dc_types, only: dp, token, string
  implicit none

  private
  public:: units
  public:: clear, deallocate, add_okay
  public:: assignment(=), operator(*), operator(/), operator(+)
  private:: units_simplify

  type units
    real(DP):: factor
    integer:: nelems
    character(TOKEN), pointer:: name(:)
    character(TOKEN):: offset
    real(DP), pointer:: power(:)
  end type units

  interface clear
    module procedure dcunitsclear
  end interface

  interface deallocate
    module procedure dcunitsdeallocate
  end interface

  interface assignment(=)
    module procedure dcunitsbuild, dcunitstostring
  end interface

  interface operator(*)
    module procedure dcunitsmul
  end interface

  interface operator(/)
    module procedure dcunitsdiv
  end interface

  interface operator(+)
    module procedure dcunitsadd
  end interface

contains

  subroutine units_simplify(u, name, power)
    type(units), intent(inout):: u
    character(*), intent(in):: name(u%nelems)
    real(DP), intent(in):: power(u%nelems)
    integer:: i, n, j, onazi
    integer:: table(u%nelems)

    if (u%nelems < 1) return
    table(:) = 0
    n = 0
    do, i = 1, u%nelems
      if (name(i) == '') cycle
      onazi = 0
      do, j = 1, i - 1
        if (name(j) == name(i)) then
          onazi = j
        endif
      enddo
      if (onazi > 0) then
        table(i) = table(onazi)
      else
        n = n + 1
        table(i) = n
      endif
    enddo
    allocate(u%name(n), u%power(n))
    u%power = 0.0_dp
    do, i = 1, u%nelems
      if (table(i) == 0) cycle
      u%name(table(i)) = name(i)
      u%power(table(i)) = u%power(table(i)) + power(i)
    enddo
    u%nelems = n
  end subroutine units_simplify

  type(units) function dcUnitsMul(u1, u2) result(result)
    type(units), intent(in):: u1, u2
    integer:: n
    character(TOKEN), allocatable:: name(:)
    real(DP), allocatable:: power(:)
    result%factor = u1%factor * u2%factor
    result%nelems = u1%nelems + u2%nelems
    result%offset = ""
    n = result%nelems
    if (n == 0) then
      nullify(result%name, result%power)
      return
    endif
    allocate(name(n), power(n))
    name = (/u1%name, u2%name/)
    power = (/u1%power, u2%power/)
    call units_simplify(result, name, power)
    deallocate(name, power)
  end function dcunitsmul

  type(units) function dcUnitsDiv(u1, u2) result(result)
    type(units), intent(in):: u1, u2
    integer:: n, n1
    character(TOKEN), allocatable:: name(:)
    real(DP), allocatable:: power(:)
    if (abs(u2%factor) < tiny(u2%factor)) then
      result%factor = sign(u1%factor, 1.0_dp) * &
        & sign(u2%factor, 1.0_dp) * &
        & huge(1.0_dp)
    else
      result%factor = u1%factor / u2%factor
    endif
    result%nelems = u1%nelems + u2%nelems
    result%offset = ""
    n = result%nelems
    if (n == 0) then
      nullify(result%name, result%power)
      return
    endif
    allocate(name(n), power(n))
    n1 = u1%nelems
    if (n1 >= 1) then
      name(1:n1) = u1%name(1:n1)
      power(1:n1) = u1%power(1:n1)
    endif
    n1 = n1 + 1
    if (n >= n1) then
      name(n1:n) = u2%name(1:u2%nelems)
      power(n1:n) = -u2%power(1:u2%nelems)
    endif
    call units_simplify(result, name, power)
    deallocate(name, power)
  end function dcunitsdiv

  type(units) function dcUnitsAdd(u1, u2) result(result)
    type(units), intent(in):: u1, u2
    type(units):: x
    result%offset = u1%offset
    result%nelems = u1%nelems
    result%factor = u1%factor + u2%factor
    x = u1 / u2
    if (x%nelems == 0) then
      nullify(result%name, result%power)   
      return
    endif
    if (all(abs(x%power(1:result%nelems)) < tiny(0.0_dp))) then
      allocate(result%name(result%nelems), result%power(result%nelems))
      result%name = u1%name
      result%power = u1%power
      return
    endif
    result%factor = 0.0
    result%nelems = -1
    result%offset = "MISMATCH"
    nullify(result%name, result%power)   
  end function dcunitsadd

  logical function add_okay(u1, u2) result(result)
    type(units), intent(in):: u1, u2
    type(units):: x
    character(STRING):: debug
    call clear(x)
    x = u1 / u2
    debug = u1
    debug = u2
    debug = x
    if (x%nelems == 0) then
      result = .true.
    else if (all(abs(x%power(1:x%nelems)) < tiny(0.0_dp))) then
      result = .true.
    else
      result = .false.
    endif
    call deallocate(x)
  end function add_okay

  subroutine dcunitsclear(u)
    type(units), intent(inout):: u
    nullify(u%name)
    nullify(u%power)
    u%factor = 1.0_dp
    u%offset = ""
    u%nelems = 0
  end subroutine dcunitsclear

  subroutine dcunitsdeallocate(u)
    type(units), intent(inout):: u
    if (associated(u%name)) deallocate(u%name)
    if (associated(u%power)) deallocate(u%power)
    u%factor = 1.0_dp
    u%offset = ""
    u%nelems = 0
  end subroutine dcunitsdeallocate

  subroutine dcunitstostring(string, u)
    !
    ! UNITS 型変数から文字型変数を生成します。
    !
    character(*), intent(out):: string
    type(units), intent(in):: u
    integer:: i, ip, npower
    character(TOKEN):: buffer
    character:: mul = '.'
    real(DP), parameter:: allowed = epsilon(1.0_dp) * 16.0

    if (u%nelems < 0) then
      string = 'error from ' // u%offset
      return 
    endif

    write(buffer, "(1pg20.12)") u%factor
    string = buffer
    if (u%nelems < 1) return

    if (abs(u%factor - 1.0) < allowed) then
      string = ""
    else if (abs(u%factor + 1.0) < allowed) then
      string = "-"
    endif

    ip = len_trim(string) + 1
    do, i = 1, u%nelems
      npower = nint(u%power(i))
      if (abs(1.0 - u%power(i)) < allowed) then
        buffer = " "
      else if (abs(npower - u%power(i)) < allowed) then
        write(buffer, "(i10)") npower
        buffer = adjustl(buffer)
      else
        write(buffer, "(1pg10.3)") u%power(i)
        buffer = adjustl(buffer)
      endif
      if (buffer == '0') cycle
      string = trim(string) // mul // trim(u%name(i)) // trim(buffer)
    enddo
    if (ip <= len(string)) string(ip:ip) = ' '
    if (string(1:1) == " ") string = adjustl(string)
    if (u%offset /= "") then
      string = trim(string) // '@' // trim(u%offset)
    endif
  end subroutine dcunitstostring

  subroutine dcunitsbuild(u, cunits)
    !
    ! 文字型変数から UNITS 型変数を生成します (constructor)。
    !
    use dcunits_com
    type(units), intent(out):: u
    character(*), intent(in):: cunits

    ! 構築中の情報、乗算対象の列として保持する。
    ! これは shift オペレータ付き単位を乗算できないことを示す。
    type elem_units
      character(TOKEN):: name
      real(DP):: power, factor
    end type elem_units
    type(elem_units), target:: ustack(100)
    integer:: ui = 1

    ! 構文単位が占める乗算対象の stack における最初の添字
    type paren_t
      real(DP):: factor
      integer:: factor_exp
      logical:: factor_inv
      integer:: power_exp
      integer:: paren_exp
    end type paren_t
    type(paren_t):: pstack(50)
    integer:: pi = 1

    ! パーサの状態遷移
    integer, parameter:: Y_INIT = 1, y_number = 2, y_name = 3, &
      & y_nx = 4, y_ni = 5, y_mul = 6, y_shift = 7
    integer:: yparse_status = y_init

    ! 字句
    integer:: ltype
    integer:: ivalue(5)
    real(DP):: dvalue
    character(TOKEN):: cvalue
    ! その他
    integer:: i
        
    ! --- 実行部 ---
    ! 初期化
    if (associated(u%name)) deallocate(u%name)
    if (associated(u%power)) deallocate(u%power)
    u%nelems = 0
    u%offset = ""
    u%factor = 1.0_dp
    if (cunits == "") return
    call dcunitssetline(cunits)
    call ustack_clear
    call pstack_clear
    yparse_status = y_init

    do
      call dcunitsgettoken(ltype, ivalue, dvalue, cvalue)
      select case(ltype)
      case (s_integer)
        select case(yparse_status)
        case (y_init, y_mul)
          pstack(pi)%factor = pstack(pi)%factor * ivalue(1)
          yparse_status = y_number
        case (y_name, y_nx)
          i = pstack(pi)%power_exp
          ustack(i:ui)%power = ustack(i:ui)%power * ivalue(1) 
          call power_next
          yparse_status = y_ni
        case (y_shift)
          u%offset = cvalue
        case default
          call error
        end select
      case (s_real)
        select case(yparse_status)
        case (y_init, y_mul)
          pstack(pi)%factor = pstack(pi)%factor * dvalue
          yparse_status = y_number
        case (y_name, y_nx)
          i = pstack(pi)%power_exp
          ustack(i:ui)%power = ustack(i:ui)%power * dvalue
          call power_next
          yparse_status = y_ni
        case (y_shift)
          u%offset = cvalue
        case default
          call error
        end select
      case (s_text)
        select case(yparse_status)
        case (y_init, y_number, y_mul)
          ustack(ui)%name = cvalue
          yparse_status = y_name
        case (y_name, y_ni)
          call ustack_grow
          call power_next
          ustack(ui)%name = cvalue
          yparse_status = y_name
        case (y_shift)
          u%offset = cvalue
        case default
          call error
        end select
      case (s_exponent)
        select case(yparse_status)
        case (y_name)
          yparse_status = y_nx
        case default
          call error
        end select
      case (s_multiply)
        select case(yparse_status)
        case (y_number, y_name)
          call factor_next
          yparse_status = y_mul
        case default
          call error
        end select
      case (s_divide)
        select case(yparse_status)
        case (y_number, y_name)
          call factor_next
          pstack(pi)%factor_inv = .true.
          yparse_status = y_mul
        case default
          call error
        end select
      case (s_shift)
        if (yparse_status == y_nx) call cancel_exp
        call units_finalize
        yparse_status = y_shift
      case (s_openpar)
        if (yparse_status == y_nx) call cancel_exp
        call pstack_push
      case (s_closepar)
        call units_finalize
        call pstack_pop
      case (s_eof)
        exit
      case default
        call error
      end select
    enddo
    if (yparse_status == y_nx) call cancel_exp
    call units_finalize

    u%nelems = ui
    u%factor = product(ustack(1:ui)%factor)
    call units_simplify(u, ustack(1:ui)%name, ustack(1:ui)%power)

  contains

    subroutine cancel_exp
      print *, operator"DCUnitsBuild: syntax error, (**) ignored"
    end subroutine cancel_exp

    subroutine error
      print *, "DCUnitsBuild: unexpected token <", &
        & trim(cvalue), "> ignored"
    end subroutine error

    subroutine power_next
      ! power_exp の終了処理
      call ustack_grow
      pstack(pi)%power_exp = ui
    end subroutine power_next

    subroutine factor_next
      ! factor_exp の終了処理
      real(DP):: factor
      i = pstack(pi)%factor_exp
      factor = product(ustack(i:ui)%factor) * pstack(pi)%factor
      if (pstack(pi)%factor_inv) then
        ustack(i:ui)%power = -ustack(i:ui)%power
        factor = 1.0_dp / factor
      endif
      ustack(i)%factor = factor
      ustack(i+1:ui)%factor = 1.0_dp
      call power_next
      pstack(pi)%factor = 1.0_dp
      pstack(pi)%factor_exp = ui
    end subroutine factor_next

    subroutine units_finalize
      call factor_next
    end subroutine units_finalize

    subroutine ustack_clear
      ui = 0
      call ustack_grow
    end subroutine ustack_clear

    subroutine ustack_grow
      if (ui >= size(ustack)) stop 'DCUnitsBuild: too many elements'
      ui = ui + 1
      ustack(ui)%name = ""
      ustack(ui)%factor = 1.0_dp
      ustack(ui)%power = 1.0_dp
    end subroutine ustack_grow

    subroutine pstack_clear
      pi = 0
      call pstack_push
    end subroutine pstack_clear

    subroutine pstack_push
      if (pi >= size(pstack)) stop 'DCUnitsBuild: too many parens'
      pi = pi + 1
      call ustack_grow
      pstack(pi)%factor_exp = ui
      pstack(pi)%factor = 1.0_dp
      pstack(pi)%factor_inv = .false.
      pstack(pi)%power_exp = ui
      pstack(pi)%paren_exp = ui
    end subroutine pstack_push

    subroutine pstack_pop
      ! factor_exp の終了処理
      call power_next
      pi = pi - 1
    end subroutine pstack_pop

  end subroutine dcunitsbuild

end module dc_units