Modern Fortran Explained Appendix E が intel fortran v14.0 で動かない件 その他

Modern Fortran Explained Appendix E のサンプルが intel fortran v14.0 で動かない件

M.Metcalf, J.Reid and M.Cohen の “Modern Fortran explained” の errata　が夏に出ました。ニュースグループの当該スレッドによると、付録Eのサンプルプログラムの訂正がもっとも重要とのことだったのですが、あいにく最新版の Intel Fortran でも実行時に異常終了してしまいます。 プログラムは Fortran2003 で導入されたオブジェクト指向の命令を用いたサンプルで、いわゆる円環状のリスト構造とその要素を挿入・削除するようなメソッドを実現したものになっています。

訂正済みのプログラム oo.f90 が下記サイトにも上げられていて、夏以降もいじった形跡があるのですが、実行時に異常終了する状況は変わっておりません。

https://groups.google.com/forum/?hl=en&fromgroups#!topic/comp.lang.fortran/aRz3HMpblTs

ftp://ftp.numerical.rl.ac.uk/pub/MRandC/edits.pdf

ftp://ftp.numerical.rl.ac.uk/pub/MRandC/oo.f90

仕方がないので、つらつらソースを眺めたところ、subroutine delete にまだ誤りがあり、call remove(item) は call remove(temp) ではないかと思いました。コメントとはやや異なり、item の指す先の中身をいじりつつ、item の指す番地は残しておくのが正しい気がします。こうすると正常終了します。またその場合、仮引数 item は intent(in) 属性をつけた方が良い事になります。 この辺は OOP 機能と挙動が変更されたポインタ引数が絡むので、これでいいのか良くわからんのですが…間違ってたらゴメンｗ

  !
  ! Delete an item: removes it from the list and deallocates it.
  !
  subroutine delete(item)
    class(anyitem), target  :: item ! <---- class(anyitem), intent(in), target  :: item
    class(anyitem), pointer :: temp
    temp => item
    call remove(item) ! <-----    call remove(temp)
    deallocate (temp)
  end subroutine

なお、コンパイル時に Fortran2003 解釈オプションをオンにする必要があります。

【悲報】 forall に死刑宣告！

全く期待はずれでダメダメちゃんな forall 命令に対して、イギリスの ISO 委員から廃止提案がなされていましたが、その後の ISO 小委員会でこの提案が承認された模様です。

ダメダメな forall ちゃん伝説に関しては以前の記事参照のこと。 http://fortran66.hatenablog.com/entry/20120629/1340905690

新しいお友達 do concurrent ちゃんが出来たので forall ちゃんはもう要らない子になったようです。

イギリスの提案 ftp://ftp.nag.co.uk/sc22wg5/N1951-N2000/N1975.txt 2f. The FORALL construct has failed in its purpose of enabling better performance and is redundant with the DO construct, especially the CONCURRENT form.

その詳しい理由 http://j3-fortran.org/doc/year/13/13-323.txt The FORALL construct and statement were added to the language in the expectation that they would enable highly efficient execution, especially on parallel processors. However, the experience with them indicates that they are too complex and have too many restrictions for compilers to take advantage of them. They are redundant with the DO CONCURRENT loop, and may of the manipulations for which they might be used may be done more efficiently by use of pointers, especially using pointer rank remapping."

投票結果 ftp://ftp.nag.co.uk/sc22wg5/N1951-N2000/N1977.txt

The subgroup had recommended approval of all of UK-10 (new
obsolesences and deletions) other than subitem 2c (obsolescing
implicit interfaces).  There were objections to subitems 1a, 1b, 2f
and 3.  After further discussion parts were voted on separately: the
set 2a, 2b, 2d and 2e was approved by 11-0-1, 1a and 1b (new
deletions) by 10-1-1 and 2f (obsolescing the FORALL construct) by
11-0-1.  A decision on subitem 3 (unobsolescing the CHARACTER* form)
was postponed.

馬鹿な子ほどかわいいと申しますが、forall ちゃんにもループ的な代入構造が１行で書けるという利点があるんです！

モボ・モガなら Eratosthenes の篩は実質３行で

Modern Fortran の動的プログラム仕様を持ちいると、エラトステネスの篩も実質３行というか本質的には１行で書けます。forallのおかげです。まぁ動的プログラムはすぐ stack overflow が出るんですけどね。

    program eratos0
      implicit none
      integer, parameter :: n = 10**4, m = int(sqrt(real(n)))
      integer ::i
      integer, allocatable ::prime(:)
      prime = [0, (i, i = 2, n)]
      forall(i = 2:m) prime(i**2::i) = 0  
      prime = pack(prime, prime /= 0)
      print *, prime
      print *, size(prime)  
      stop
    end program eratos0

実行結果

Steve Lionel による forall と do concurrent に関するまとめを含むスライド

http://software.intel.com/sites/default/files/comment/1717617/parallel-programming-features.pdf

上記PDF中で、do concurrent には mask 機能が無いと書いてありますが、"Modern Fortran explained" を見る限りはあるようです。また Intel Fortran でも利用できます。ループ制御変数のマスク機能が forall 最適化を邪魔する元凶だったので、何らかの制限がついたのかもしれません。詳細不明です。

Intel Fortran ver.14.0.1.139 出る！

最近、Intel Fortran の Update が出ました。コンパイラの方はあまり変化が無いようです。バグ修正くらいでしょうか？

メモ帳： Lapack95 ライブラリは手作業で指定！

Intel Visual Fortran コンパイラには LAPACK 利用のオプションがありますが、これはオリジナルの FORTRAN77 版には対応していますが、Fortran95 のラッパーの方には対応していません。

適切なディレクトリとファイルを手動で指定する必要があります。このほかに通常の LAPACK 使用オプションも指定。

メモ帳： ユーザー定義派生型変数入出力

Intel Fortran v.14 から対応になった non-default derived-type input/output の練習。Modern Fortran Explained 17.2 のサンプルを動くように適宜補ったもの。

    module m_sub
      implicit none
      type :: t_person
        character (len = 20) :: name
        integer :: iage
      contains
        procedure :: pwf
        generic   :: write(formatted) => pwf
      end type t_person    
    contains  
      subroutine pwf(dtv, unit, iotype, vlist, iostat, iomsg)
        class (t_person)   , intent(in) :: dtv
        integer            , intent(in) :: unit 
        character (len = *), intent(in) :: iotype 
        integer            , intent(in) :: vlist(:)
        integer            , intent(out) :: iostat
        character (len = *), intent(inout) :: iomsg
        !
        character (len = 9) :: pfmt
        write(pfmt, '(a, i0, a, i0, a)') '(a', vlist(1), ',i', vlist(2), ')'   
        write(unit, pfmt, iostat = iostat) dtv%name, dtv%iage
        return
      end subroutine pwf  
    end module m_sub
    
    program DT
      use m_sub 
      implicit none
      integer :: id, members
      type (t_person) :: chairman
      id      = 12345 
      members = 666
      chairman = t_person('Thatcher', 50)
      write(*, '(i8, 2x, dt(10, 6), i5)') id, chairman, members
    end program DT

実行結果

メモ帳： Merge Sort

安定ソートが欲しくなったので勉強。マージするところが間接的に再帰になっていますが、これは再帰指定しなくてもよいものなのか否や？

私は FORTRAN77 でも安心のシェルソート大好き人間です。ギャップ関数さえ選べば速さは十分なんですが、安定ソートじゃないんですよね。

シェルソートのギャップ関数に関しての以前の記事をご覧ください。http://fortran66.hatenablog.com/entry/20100226/1267120842

    module m_msort
      implicit none
    contains
      recursive function msort(x) result(res)
        real, allocatable :: res(:)
        real, intent(in)  :: x(:)
        integer :: k
        if (size(x) < 2) then
          res = x
        else
          k = size(x) / 2
          res = merge_dat( msort(x(:k)), msort(x(k + 1:)) )
        end if    
        return 
      end function msort
      
      function merge_dat(x1, x2) ! indirectly recursive ...
        real, allocatable :: merge_dat(:)
        real, intent(in) :: x1(:), x2(:)
        integer :: i, k1, k2
        allocate( merge_dat(size(x1) + size(x2)) )
        i  = 1
        k1 = 1
        k2 = 1
        do  
          if (k1 > size(x1)) then
            merge_dat(i:) = x2(k2:)
            exit 
          end if  
          if (k2 > size(x2)) then
            merge_dat(i:) = x1(k1:)
            exit  
          end if
          if (x1(k1) > x2(k2)) then 
            merge_dat(i) = x2(k2)
            k2 = k2 + 1
          else
            merge_dat(i) = x1(k1)
            k1 = k1 + 1
          end if    
          i = i + 1
        end do
        return
      end function merge_dat  
      
    end module m_msort
    
    program merge
      use m_msort
      implicit none
      real, allocatable :: x(:)
      allocate( x(10**2) )
      call random_number(x)
      x = msort(x)
      print *, x
      print *, all( x <= eoshift(x, 1, huge(x)) ) ! check
      stop 
    end program merge

実行結果