fortran66のブログ

fortran について書きます。

【メモ帳】エラトステネスの篩

Sieve of Eratosthenes

julia と fortran でおおむね同じなプログラムで比較。10**9 以下の素数をエラトステネスの篩で求める。

Julia 1.5.1

julia さん速くて Fortran -O3 と揺らぎの範囲で大差ない。

@time pr = prim(10^9)
33.018580 seconds (5 allocations: 14.901 GiB, 0.30% gc time)
function prim(n::Int64)
    tab = Vector(1:n)
    tab[1] = 0
    for i = 2:floor(Int64, sqrt(n))
        if tab[i] != 0
            for j = i*i:i:n
                tab[j] = 0
            end
        end
    end
    filter(!iszero, tab)
end
50847534-element Array{Int64,1}:
         2
         3
         5
         7
        11
        13
        17
        19
        23
        29
        31
        37
        41
         ⋮
 999999667
 999999677
 999999733
 999999739
 999999751
 999999757
 999999761
 999999797
 999999883
 999999893
 999999929
 999999937

Intel Fortran 19.1.2

/heap-arrays100000 として大きい配列を stack から heap に移さないと stack overflow する。

int64     50847534   33.68400
int32     50847534   27.19900

今の時代の計算機はメモリー移動が律速段階だから 32bit にすると速い。もしくは Fortran の default 長が 32bit のためかもしれない。64 bit 整数対応は後付け。

module prim_m
    use, intrinsic :: iso_fortran_env
    implicit none
contains
    function prime(n) result(ires)
        integer, intent(in) :: n
        integer(int64), allocatable :: ires(:)
        integer(int64) :: i, itab(n)

        forall (i = 1:n) itab(i) = i
        itab(1) = 0
        do i = 2, int(sqrt(real(n)))
            if (itab(i) /= 0) itab(i*i::i) = 0
        end do
        ires = pack(itab, itab /= 0)
    end function prime
end module prim_m


program prim
    use prim_m
    implicit none
    integer(int64), allocatable :: ip(:)
    integer(int64) :: ic0, ic1, irate
    call system_clock(ic0)
    ip = prime(10**9)
    call system_clock(ic1, irate)
    print *, size(ip), (ic1 - ic0) / real(irate)
end program prim

1から始める Juliaプログラミング

1から始める Juliaプログラミング

数値計算のためのFortran90/95プログラミング入門(第2版)

数値計算のためのFortran90/95プログラミング入門(第2版)

  • 作者:牛島 省
  • 発売日: 2020/01/28
  • メディア: 単行本(ソフトカバー)

Fortran90/95プログラミング

Fortran90/95プログラミング

Fortran ハンドブック

Fortran ハンドブック

【訃報】Walter Brainerd 氏死去

おくやみ

六月に亡くなられていたようです。

www.reddit.com

www.fortran.com

Guide to Fortran 2008 Programming

Guide to Fortran 2008 Programming

Programmer's Guide to Fortran 90

Programmer's Guide to Fortran 90

Programmer's Guide to F

Programmer's Guide to F

Fortran 77: Fundamentals and Style

Fortran 77: Fundamentals and Style

【メモ帳】Arjen Markus ゲスト出演

Fortran Pointer Vs Allocatable With Arjen Markus

Fortran系Youtuber? な Everything Functional さんの動画に Arjen Markus さんがゲスト出演。動的メモリー確保用の allocatable と pointer について語っています。


Fortran Pointer Vs Allocatable With Arjen Markus

冒頭紹介でアルジェン・マーカス(英語読み?)の個人情報が開示されているのでストーカーの皆さんは( ..)φメモメモ

参考

allocatable と pointer に関しては、以下の John Reid の講演も参考になります。

fortran66.hatenablog.com

【メモ帳】画像変換

屈辱の画像変換w

Fortran で jpg や png を読み書きするのに、Python スクリプトを内部生成してテンポラリな bmp ファイルを経由して I/O するという屈辱的な方法を使ってみます。まぁどうせどんなファイルも読み込めば BMP になるんだし~手段をえらばず。

python script を生成して実行。bmp ファイルに/から変換する

        subroutine conv_bmp(fn_in, fn_out)
            character(*), intent(in) :: fn_in, fn_out
            character(*), parameter :: fn_tmp = 'my_python_tmp.py'
            integer :: iw
            open(newunit = iw, file = fn_tmp)
            write(iw, '(a)') 'from PIL import Image'
            write(iw, '(*(a))') 'img = Image.open("', trim(fn_in), '").convert("RGB")'
            write(iw, '(*(a))') 'img.save("', trim(fn_out), '")'
            close(iw)
            call execute_command_line('python3 '//fn_tmp)
            call delete_file(fn_tmp)
        end subroutine conv_bmp

実行例

 image size         535         210

入力 1966.png

f:id:fortran66:20201002030034p:plain

出力 1966_without_red.jpg

赤色を消します。符号なし 8 bit 整数の代用として文字を使っています。

bmp%rgb(:,:)%r = achar(0)

f:id:fortran66:20201002030049j:plain

open(iw, fn)
close(iw, status = 'delete')

でテンポラリファイルを消しています。

write(iw, '(*(a))') 'img = Image.open("', trim(fn_in), '").convert("RGB")'

Python の pillow を使っています。Fortran 側の BMP I/O が 24 bit RGB 専用なので、読み込む時に .convert("RGB") をつけておく必要があるようです。

    module m_bmp
        use, intrinsic :: iso_fortran_env
        implicit none
        type :: t_bmp_file_header
            sequence  
            integer(int16) :: bfType = transfer('BM', 0_int16) ! BitMap
            integer(int32) :: bfSize          ! file size in bytes
            integer(int16) :: bfReserved1 = 0 ! always 0
            integer(int16) :: bfReserved2 = 0 ! always 0
            integer(int32) :: bfOffBits
        end type t_bmp_file_header
        !
        type :: t_bmp_info_header
            sequence
            integer(int32) :: biSize     = Z'28' ! size of bmp_info_header ; 40bytes 
            integer(int32) :: biWidth
            integer(int32) :: biHeight
            integer(int16) :: biPlanes   = 1 ! always 1
            integer(int16) :: biBitCount
            integer(int32) :: biCompression = 0 ! 0:nocompression, 1:8bitRLE, 2:4bitRLE, 3:bitfield
            integer(int32) :: biSizeImage
            integer(int32) :: biXPelsPerMeter = 3780 ! 96dpi
            integer(int32) :: biYPelsPerMeter = 3780 ! 96dpi 
            integer(int32) :: biClrUsed       = 0
            integer(int32) :: biClrImportant  = 0 
        end type t_bmp_info_header
      !
        type :: t_rgb
            sequence
            character :: b, g, r
        end type t_rgb  
      !
        type :: t_bmp
            type(t_bmp_file_header) :: file_header
            type(t_bmp_info_header) :: info_header
            type(t_rgb), allocatable :: rgb(:, :)
        contains
            procedure :: rd => rd_py      
            procedure :: wr => wr_py
            procedure :: rd_bmp      
            procedure :: wr_bmp
        end type t_bmp    
    contains   
        subroutine rd_py(bmp, fn)
            class(t_bmp), intent(out) :: bmp
            character(len = *), intent(in) :: fn
            character(*), parameter :: fn_tmp = 'my_bmp_rd_tmp'
            call conv_bmp(fn, fn_tmp//'.bmp')
            call bmp%rd_bmp(fn_tmp)
            call delete_file(fn_tmp//'.bmp')
        end subroutine rd_py
        
        
        subroutine wr_py(bmp, fn)
            class(t_bmp), intent(in out) :: bmp
            character(len = *), intent(in) :: fn
            character(*), parameter :: fn_tmp = 'my_bmp_wr_tmp'
            call bmp%wr_bmp(fn_tmp)
            call conv_bmp(fn_tmp//'.bmp', fn)
            call delete_file(fn_tmp//'.bmp')
        end subroutine wr_py
        
        
        subroutine rd_bmp(bmp, fn)
            class(t_bmp), intent(out) :: bmp
            character(len = *), intent(in) :: fn
            type(t_bmp_file_header) :: bmp_file_header 
            type(t_bmp_info_header) :: bmp_info_header
            integer :: ix, iy, nx, ny, ir
            character :: dummy
            open(newunit = ir, file = fn//'.bmp', access = 'stream', status = 'unknown')
            read(ir) bmp_file_header, bmp_info_header
            nx = bmp_info_header%biWidth
            ny = bmp_info_header%biHeight
            allocate( bmp%rgb(nx, ny) )
            read(ir) (bmp%rgb(:, iy), (dummy, ix = 1, mod(nx, 4)), iy = 1, ny)
            close(ir)
        end subroutine rd_bmp  
        
    
        subroutine wr_bmp(bmp, fn)
            class(t_bmp), intent(in out) :: bmp
            character(len = *) :: fn
            integer :: ix, iy, nx, ny, iw
            nx = size(bmp%rgb, 1)
            ny = size(bmp%rgb, 2)
            bmp%file_header%bfSize      = 14 + 40 + 0 + (3 * nx + mod(nx, 4)) * ny
            bmp%file_header%bfOffBits   = 14 + 40
            bmp%info_header%biWidth     = nx
            bmp%info_header%biHeight    = ny
            bmp%info_header%biBitCount  = 24 
            bmp%info_header%biSizeImage = (3 * nx + mod(nx, 4)) * ny
            open(newunit = iw, file = trim(fn)//'.bmp', access = 'stream', status = 'unknown')
            write(iw) bmp%file_header, bmp%info_header
            write(iw) (bmp%rgb(:, iy), (achar(0), ix = 1, mod(nx, 4)), iy = 1, ny)
            close(iw)
        end subroutine wr_bmp  

        
        subroutine conv_bmp(fn_in, fn_out)
            character(*), intent(in) :: fn_in, fn_out
            character(*), parameter :: fn_tmp = 'my_python_tmp.py'
            integer :: iw
            open(newunit = iw, file = fn_tmp)
            write(iw, '(a)') 'from PIL import Image'
            write(iw, '(*(a))') 'img = Image.open("', trim(fn_in), '").convert("RGB")'
            write(iw, '(*(a))') 'img.save("', trim(fn_out), '")'
            close(iw)
            call execute_command_line('python3 '//fn_tmp)
            call delete_file(fn_tmp)
        end subroutine conv_bmp
        

        subroutine delete_file(fn)
            character(*), intent(in) :: fn
            integer :: iw
            open(newunit = iw, file = fn)
            close(iw, status = 'delete')
        end subroutine delete_file
    end module m_bmp
    

    program imgage
        use, intrinsic :: iso_fortran_env
        use :: m_bmp
        implicit none
        type(t_bmp) :: bmp
        call bmp%rd('1966.png')
        print *, 'image size', shape(bmp%rgb)
        bmp%rgb(:, :)%r = achar(0)
        call bmp%wr('1966_without_red.jpg')
    end program imgage

R2-10-2 微修正 delete_file 追加 R2-10-3 微修正 tmp file 名 parameter 化

生成 Deep Learning ―絵を描き、物語や音楽を作り、ゲームをプレイする

生成 Deep Learning ―絵を描き、物語や音楽を作り、ゲームをプレイする

  • 作者:David Foster
  • 発売日: 2020/10/05
  • メディア: 単行本(ソフトカバー)

最近オライリー本が白くなってマイケル・ジャクソンみたい。中の紙も漂白紙になって目が痛い。ホワイトウォッシュだなw

【乞食速報】Humble Bundle No Starch Press その他寝言

HUMBLE BOOK BUNDLE: LEARN TO CODE THE FUN WAY BY NO STARCH PRESS

www.humblebundle.com

社会主義リアリズム芸術としての杜甫

杜甫の詩は社会批判になると社会主義リアリズム的でイライラします。

各種障害

Google, Microsoft, Apple と日替わりで連日障害を起こしていますが、とある外国の仕業である可能性が無くも無いアルヨ!ポコペン

杜甫―偉大なる憂鬱

杜甫―偉大なる憂鬱

「ズベズダ(Zvezda)」で空気漏れ

www.afpbb.com

www.sekaiseifuku-zzz.com

【寝言】ニュース

NEC VECTOR

www.nextplatform.com

Fujitsu A64FX

ノードあたりのメモリーが少ない。 それが問題にならなければ x86 と互角以上?

EU の夢

イギリスは当然ながらおらほの ARM だろうが、大陸EUが RISC-V に走っていて草。

FORTRAN vs PASCAL

https://pbs.twimg.com/media/EjBjsL3WsAAYgsO?format=jpg&name=small https://pbs.twimg.com/media/EjBjtAtWsAEDc1K?format=jpg&name=small

【メモ帳】ランド研究所の乱数

Rand's random numbers

WSJランド研究所の百万桁乱数本に関する記事が出ています。以前、ハーマン・カーンからみで出てきた乱数表だと思いますが、電流の雑音を数値化してパンチカードに打ち出して乱数に直したようですが、検証した所一部に異常値があるので、パンチカードの順序が入れ替わったにではと推測しているようです。

www.wsj.com

ネタ元:

コメント欄によると、無料で件の百万桁乱数本が落とせるようです。 www.rand.org

fortran66.hatenablog.com

Mathematical Methods for Digital Computers: v. 1

Mathematical Methods for Digital Computers: v. 1

  • 作者:RALSTON MATHEMA
  • 発売日: 1960/12/01
  • メディア: ハードカバー