fortran66のブログ

fortran について書きます。

【メモ帳】WebGPU とやら

M1 mac WebGPU

斯様な記事が出ていたので、試してみました。結果表示を信じるなら確かに約900GFlopsを出しています。

gigazine.net

WebGPUもよく分かりませんが、JavaScriptもよく分かりません。Flops導出式が正しいのかもよく分かりません。

出力結果

best: 933.69 gflops
[numthreads(4, 16, 1)]
compute void main(constant float4[] A : register(u0),
                  constant float4[] B : register(u1),
                  device float4[] C : register(u2),
                  float3 threadID : SV_DispatchThreadID) {
  uint m = uint(threadID.x);
  uint n = uint(threadID.y);

以下略

しかし秘密主義のアップルのせいでよく分からないGPUも、アップル絡みのWebGPUに関してはブラウザ上のJavaScriptでハードウェアに近いAPIを叩けるという奇妙さ。

GPUの他に、ニューラルエンジンというのもあるけれども、これも何が何やらよく分かりません。マニュアルを見ていると整数1024byteのSIMD演算がAPIがあるように見受けられるけれど、よく分かりません。SwiftかObjective-Cしか例題もなし。CのAPIがあればfortranからも呼べようけども・・・

【寝言】Hey! Siri!

中国がコロナ肛門検査

www.jiji.com

中国駐在の米国務省職員が「誤って」検査を受けたと報じた。

米外交官にコロナ肛門検査? 中国が報道否定
2021年02月26日11時32分

米外交官にコロナ肛門検査? 中国が報道否定
【北京AFP=時事】中国政府は25日、中国で米国務省職員が新型コロナウイルス検査のため肛門からの検体採取を受けたとした米メディアの報道を否定し、米外交官にこうした検査を求めたことは一度もないと説明した。

 国内での感染拡大がほぼ抑えられている中国は先月、ウイルスは消化器官内でより長く残存するため、検体を採取する部位は通常の喉や鼻ではなく肛門の方が効果的である可能性を指摘していた。

 米オンラインメディアのバイスと米紙ワシントン・ポストは米当局者の話として、外交官が肛門検査の対象外とされているのにもかかわらず、中国駐在の米国務省職員が「誤って」検査を受けたと報じた。

 中国側はこの報道を否定。同国外務省の趙立堅報道官は25日の定例記者会見で、「中国が国内駐在の米外交職員に肛門からの検体採取を受けるよう要求したことは一度もない」と断言した。

 中国当局は、感染リスクが高いとされる人々に対し、肛門からの検体採取を通じた検査を実施してきた。こうした人々には、感染が起きた地域の住民のほか、国外からの渡航者の一部も含まれる。1月には、北京市で小規模の集団感染が発生した際にこの検査が行われたとの国営メディアの報道がソーシャルメディア上で反響を呼び、中国版ツイッターの「微博(ウェイボー)」では笑いと恐怖が入り交じった反応が相次いだ。【翻訳編集AFPBBNews】
〔AFP=時事〕

f:id:fortran66:20210228163840j:plain

Hey! Siri! は勘弁してほしい。数値解析(Numerical Analysis) をヌメリ・アナルと呼ぶ趣き(全然違う!)。

そもそも siri は植木等ばりにおよびでない時に出てきて、頓珍漢な答えしか返さないので、マイクロソフトのイルカ並みに殺意が沸く。

Introduction to Numerical Analysis (Texts in Applied Mathematics, 12)

Introduction to Numerical Analysis (Texts in Applied Mathematics, 12)

  • 作者:Stoer, J.
  • 発売日: 2010/02/19
  • メディア: ペーパーバック

【メモ帳】メモ

M1 Mac メモ

intel devcloud for oneAPI と visual studio code

Intel devcloud に ARM 版 visual studio code で接続して、Modern Fortran Extensions を入れようとすると、間接的に利用している c++ extensions が途中で引っかかってサーバーから切れてしまう。 intelvisual studio code の場合は、切れても何回か reload すれば突破できるが、ARM 版の場合 reload してもうまくゆかない。

Modern Fortran Extensions は、ローカルでは問題ないように見えています。

よく分からないし、あまり調べていない。

どうも昨夜はサーバーの調子が悪かっただけのようで、ちゃんと指示通りに port forwarding してつなげばとりあえずは使えました。ただ左側のメニューアイコンに intel oneAPI の記号が現れないし時々ゆえ無くして切れるので、何かは問題な模様。 手動で sample browser for oneAPI を入れればよい模様です。

devcloud.intel.com

電池の減り方

タブレットの電池の減り方に似ている。丸一日触らずスリープ状態で放置していても、電池は減らない。ブラウザでネットを見ているとゆるゆる減ってゆく。Bluetooth を使うとグングン減る。

ふたを開けるともうロック画面が出ていて、指紋認証スイッチですぐロック解除されるところはタブレットの様。

parallel desktop 上の Windows 10 に visual studiointel fortran をインストールして実行

parallel desktop と ARM Windows10 のインストールはこちらのページを参考にしました。ありがたいです。

pochilog.jp

WSL は WSL1 も WSL2 もうまく動かなかったので、(まだ parallel desktop で virtual machine の仮想化の多重化が指定できないようになっている)、windows 側に intel fortran をインストールしました。先に visual studio 2019 を C++ desktop SDK と一緒にインストールしておく必要があります。oneAPI の HPC toolkit から MPI と Fortran のみをインストールしました。特に問題なく進みます。

コンパイルできるのは 32bit 版(x86)のプログラムだけで、64bit 版(x64) の方はエラーが出てコンパイルできませんでした。x86 版は少し前から coarray 対応を削除されているので、こちらでも coarray は試せません。

mac os 側で試したのと同じプログラムでベンチマークしてみますと、Rosetta 2 の倍くらいの時間がかかって gfortran と同程度の実行時間になっています。なお intel cpu では、一般的に x86 版は x64 版に比べてちょこっと遅い事が多いです。

2021-02-28  01:22:25.4   cpu_time =       0.00    (sec) ::make matrix
2021-02-28  01:22:25.5   cpu_time =      0.625E-01(sec) ::solve Ax = b
 difference (rms) =  5.967297611830250E-007
2021-02-28  01:22:27.6   cpu_time =       2.20    (sec) ::normal end
Press any key to continue . . .

f:id:fortran66:20210228013029p:plain

Apple IIZ80 ボードを挿して CP/M を動かすようでいいですね。

【メモ帳】abstract type の中に concrete な routine をぶら下げる

記事を移動し独立させました

昔のプログラムが、今の ifort を通らないので修正。

元々 abstract type に concrete なルーチンをくっつけていて違和感があったので、書き直してみました。

どう直せばいいのか少し苦労しましたが、メンバー要素として type bound procedure pointer にすることで解決しました。

fortran66.hatenablog.com

    module m_base
        implicit none
        private
        integer, parameter, public :: kd = kind(1.0d0)
!
        type, abstract, public :: t_matrix
            procedure(cg_sub), pointer :: cg => cg_sub
        contains
            procedure(prd_fun), deferred :: product_abst
            generic, public :: operator(*) => product_abst
        end type t_matrix
    
        abstract interface
            function prd_fun(mat, vec) result(res) 
                import :: t_matrix, kd 
                class(t_matrix), intent(in) :: mat
                real(kd), intent(in) :: vec(:) 
                real(kd) :: res(size(vec))
            end function prd_fun
        end interface     

    contains

        subroutine cg_sub(A, b, x) ! Conjugate Gradient method
            class(t_matrix), intent(in) :: A
            real(kd), intent(in) :: b(:)
            real(kd), intent(in out) :: x(:)
            real(kd) :: alpha, beta, r2n, r2d, eb2, p(size(b)), q(size(b)), r(size(b))
            integer :: i
            r = b - A * x
            p = r
            r2n = dot_product(r, r)
            eb2 = epsilon(b) * dot_product(b, b)
            do i = 1, size(b)
                if ( r2n < eb2 ) exit
                r2d = r2n
                q = A * p
                alpha = dot_product(p, r) / dot_product(p, q)
                x = x + alpha * p
                r = r - alpha * q
                r2n = dot_product(r, r)
                beta = r2n / r2d
                p = r + beta * p            
            end do
        end subroutine cg_sub
    end module m_base

    module m_dense
        use m_base
        implicit none
        private
        type, extends(t_matrix), public :: t_dense_matrix
            real(kd), allocatable :: x(:, :)
        contains
            procedure, public :: product_abst => product_dense  
        end type t_dense_matrix

    contains
  
        function product_dense(mat, vec) result(res)
            class(t_dense_matrix), intent(in) :: mat
            real(kd), intent(in) :: vec(:)
            real(kd) :: res(size(vec))
            res = matmul(mat%x, vec)
        end function product_dense
    end module m_dense
  
    module m_sparse
        use m_base
        implicit none
        private
        type :: t_sparse_vec
            integer , allocatable :: idx(:)
            real(kd), allocatable :: x(:)
        end type t_sparse_vec
    
        type, extends(t_matrix), public :: t_sparse_matrix
            type(t_sparse_vec), allocatable :: row(:) 
        contains
            procedure :: product_abst => product_sparse  
        end type t_sparse_matrix

    contains
    
        function product_sparse(mat, vec) result(res)
            class(t_sparse_matrix), intent(in) :: mat
            real(kd), intent(in) :: vec(:)
            real(kd) :: res(size(vec))
            real(kd) :: tmp
            integer :: i, k
            do i = 1, size(mat%row)
                tmp = 0.0_kd
                do k = 1, size(mat%row(i)%idx)
                    tmp = tmp + mat%row(i)%x(k) * vec(mat%row(i)%idx(k))
                end do    
                res(i) = tmp
            end do 
        end function product_sparse
    end module m_sparse
  
    program CG
        use m_base
        use m_dense
        use m_sparse
        implicit none
        type(t_dense_matrix ) :: a
        type(t_sparse_matrix) :: q
        real(kd), allocatable :: x(:), y(:), b(:)
        integer :: i, j, n = 1000  
    !
        allocate( a%x(n, n), b(n), x(n) ) ! dense
        allocate( q%row(n)       , y(n) ) ! sparse
    !
        do i = 1, n
            b(i) = 1.0_kd
            do j = i, n ! x = (1, 0, 0, ... , 0)t
                a%x(i, j) = real( 2 * MIN(i, j) - 1, kd )  
                a%x(j, i) = a%x(i, j)
            end do
        end do      
        do i = 1, n                             ! copy dense-mat to sparse-mat 
            q%row(i)%idx = [(j, j = 1, n)]  
            q%row(i)%x   = a%x(i, :)
        end do      
    !
        call random_seed()
        call random_number(x) ! starting vector
        y = x
        call a%cg(b, x)
        call q%cg(b, y)
        print *, 'dense  matrix |a * x - b| =', norm2(a * x - b) ! norm2 : f2008 
        print *, 'sparse matrix |q * y - b| =', norm2(q * y - b)
        print *, 'x==y?', all(x==y)
    end program CG

【メモ帳】M1 Mac で Fortran

天皇陛下のお誕生日を祝して M1 mac 購入

M1 mac を購入したので、Fortran をインストールしてみました。私は mac は生まれて初めてで右も左も分からぬ状態なので、話半分で読んでください。買ったのは一番安い macbook air 8G/256G/7GPU ので色は給食のアルマイトの食器の様なねずみ色のやつです。名前は Air ですが、ずっしり重くてどちらかというと Heavy Metal という感じです。

約1.29 kg~

最新 Apple Mac mini Apple M1 Chip (8GB RAM, 512GB SSD)

最新 Apple Mac mini Apple M1 Chip (8GB RAM, 512GB SSD)

  • 発売日: 2020/11/17
  • メディア: Personal Computers

約634g~

三種の Fortran

Arm 用 gfortran, Rosetta エミュレータintel 用 gfortran および oneAPI の intel fortran の三種が動くことを確かめました。インストールはいずれもインストーラを動かすだけでできました。

簡単なプログラムで見たところ intel fortran の場合、エミュレーションで動いているにもかかわらず、第四世代 Haswell i7 の4倍くらいのスピードが出ています。同じプログラムでみると ARM 版 gfortran は x86_64 版 gfortran の 1.5 倍くらいの速さで、intel fortran と比較するとプログラムによって速かったり遅かったりします。

ARM ネイティブ gfortran が Rosetta のエミュレーションの intel fortran より遅いのは意外ですが、ネットを流れているベンチマークを見ると、ARM 版 gcc は clang より平均2倍くらい遅いようなので、 GNU コンパイラの問題なのかもしれません。また元々 gfortran は intel fortran より遅く5割増し位な感じなので、合わせるとそんなところかもしれません。

なお、Appleintel fortran は昔から coarray をサポートせず oneAPI 版でもサポートされません。

インストール・メモ

intel fortran

Jacob Williams さんのツイートで intel fortran が動くことが分かっていました。

以下のページから HPC base kit 用のオンライン版インストーラを取ってこれます。登録が必要ですが無料です。 fortran コンパイラだけでいいなら oneAPI base kit の方は要りません。ただ Intel MKL は base kit の方に入っているので、これも入れるならば base kit 用のインストーラも必要になります。

software.intel.com

Mac 版は MPI その他が無く(多分そのせいで coarray もない... MPI 上に実装されているので)、選択できるのは C++Fortranコンパイラだけです。Fortran コンパイラだけインストールしました。

コマンドラインXcode というものが必要になるのですが、ifort でコンパイルしようとすると警告が出て自動でインストールしてくれます。

メモ コンパイラにパスを通すスクリプトを先に動かす必要があります。

 . /opt/intel/oneapi/setvars.sh 

ただライブラリへのパスが通っておらず、ネットを検索すると手動でライブラリ位置を指定せよと出てきました。Xcode の前に intel fortran をインストールしたせいかもしれません。

ifort   -L/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/lib   xxxx.f90

めんどくさいので .zshrc に alias を書いてみました。

alias ifort='ifort   -L/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/lib'

gfortran

奥村先生のページを参考にしました。 oku.edu.mie-u.ac.jp

Homebrew というドブロク造りをイメージしたパッケージマネージャーを x86_64 用と arm64 用それぞれに、ページ画面に表示されたコマンド1行をターミナルにコピペ&実行でインストールして、その後それぞれで brew install gfortran を行えば gfortran がインストールされます。

ARM 版を優先したい場合、実行パスの優先順序を工夫する必要があります。

visual studio code

insider という所に行くと ARM 版が落とせます。

code.visualstudio.com

よく分からない形式の圧縮ファイルをクリックすると、実行ファイルらしきものが出てくるので、ええいままよとそれをアプリケーションのフォルダーにコピーしてみました。これでいいのかよく分かりません。クリックすると vs code が起動します。

visual studio code 起動後は Modern FortranFortran IntelliSense の拡張機能を入れました。何か言われるままにクリックして、文法解釈用の fortls も別途インストールした気がします。漫然とクリックしたのでよく覚えていません。とりあえず問題なく動いています。

実行例

とりあえず昔のプログラムを動かしてみました。

fortran66.hatenablog.com

    module m_stamp
        ! use portlib ! MS-FPS4.0
        implicit none
        logical :: first = .true.
        real(kind(0.0d0)) :: t0, t1
    contains
        subroutine stamp(text) ! print time stamp
            character(len = *), intent(in) :: text
            character(len =  8) :: date
            character(len = 10) :: time
            if (first) then 
                first = .false.
                call cpu_time(t0) ! t0 = timef() ! MS-FPS4.0
            end if 
            call cpu_time(t1)   ! t1 = timef() ! MS-FPS4.0
            call date_and_time(date, time)
            print '(6a, 6a, a, g15.3, 2a)', &
              date(1:4), '-', date(5:6), '-', date(7:8), '  ', &
              time(1:2), ':', time(3:4), ':', time(5:8), '  ', &
              ' cpu_time =', t1 - t0, '(sec) ::', text 
        end subroutine stamp  
    end module m_stamp    
  
    module m_cg
        implicit none
        integer, parameter :: kd = kind(0.0d0) !double precision
    contains
        subroutine cg(a, b, x) ! conjugate gradient method
            real(kd), intent(in) :: a(:, :), b(:)
            real(kd), intent(in out) :: x(:)
            real(kd) :: alpha, beta, r2n, r2d, eb2, p(size(b)), q(size(b)), r(size(b))
            integer :: i
            r = b - matmul(a, x) 
            p = r
            r2n = dot_product(r, r)
            eb2 = epsilon(b) * dot_product(b, b)
            do i = 1, size(b)
                if ( r2n < eb2 ) exit
                r2d = r2n
                q = matmul(a, p)
                alpha = dot_product(p, r) / dot_product(p, q)
                x = x + alpha * p
                r = r - alpha * q
                r2n = dot_product(r, r)
                beta = r2n / r2d
                p = r + beta * p              
            end do
        end subroutine cg
    end module m_cg
  
    program cg_main
        use m_cg
        use m_stamp
        implicit none
        integer, parameter :: ns = 2000 
        real(kd) :: a(ns, ns), b(ns), x(ns)
        integer :: i, j
        call stamp('make matrix')
        do i = 1, ns
            b(i) = 1.0_kd
            do j = i, ns ! givens matrix,  x = (1, 0, 0, ... , 0)t
                a(i, j) = real( 2 * min(i, j) - 1, kd )  
                a(j, i) = a(i, j)
            end do
        end do
        call stamp('solve Ax = b')
  ! 
        call random_seed() 
        call random_number(x) ! starting vector
        call cg(a, b, x)     
        x = matmul(a, x) - b
        print *, 'difference (rms) =', sqrt( dot_product(x, x) )
        call stamp('normal end')
    end program cg_main

M1 Mac

intel fortran (Rosetta2)

gfortran (ARM64)

gfortran (Rosetta2)

[a] M1:~/fortran% ifort -O2 cg.f90                               
[a] M1:~/fortran% ./a.out         
2021-02-24  02:07:14.7   cpu_time =      0.500E-05(sec) ::make matrix
2021-02-24  02:07:14.7   cpu_time =      0.182E-01(sec) ::solve Ax = b
 difference (rms) =  4.684681968851512E-007
2021-02-24  02:07:15.6   cpu_time =      0.945    (sec) ::normal end

[a] M1:~/fortran% gfortran -O2 cg.f90                            
[a] M1:~/fortran% ./a.out            
2021-02-24  02:07:24.3   cpu_time =      0.700E-05(sec) ::make matrix
2021-02-24  02:07:24.3   cpu_time =      0.137E-01(sec) ::solve Ax = b
 difference (rms) =   6.6372870430592568E-007
2021-02-24  02:07:26.1   cpu_time =       1.80    (sec) ::normal end

[a] M1:~/fortran% arch -x86_64 /usr/local/bin/gfortran -O2 cg.f90
[a] M1:~/fortran% ./a.out                                        
2021-02-24  02:07:34.2   cpu_time =      0.120E-04(sec) ::make matrix
2021-02-24  02:07:34.2   cpu_time =      0.218E-01(sec) ::solve Ax = b
 difference (rms) =   4.3568034320423338E-007
2021-02-24  02:07:36.4   cpu_time =       2.26    (sec) ::normal end

i7 4770

Windows x64 intel fortran

2021-02-24  02:08:36.8   cpu_time =       0.00    (sec) ::make matrix
2021-02-24  02:08:36.9   cpu_time =      0.469E-01(sec) ::solve Ax = b
 difference (rms) =  4.077532388488979E-007
2021-02-24  02:08:40.8   cpu_time =       3.95    (sec) ::normal end
続行するには何かキーを押してください . . .

追記: i7 4770 の計算時間が昔に比べて長くなって計算スピードが遅くなっていますが、実は MKL 付属の linpack benchmark なども遅くなっています。多分 intel chip の specter とかの対策のせいではないかと疑っているのですが面倒なので確かめていません。気のせいかもしれません。

【寝言】細い月

二日月?

今日は暖かな日和の日本晴れで、夕焼け空に細い月が出ていました。 暗くなってくると地球照も見えて中々よい眺めでした。

f:id:fortran66:20210215003349j:plain

野鳥のラップバトル

最近、小鳥の縄張り争いが繰り広げられていて、一羽がひとしきりギーギーと耳障りな声で啼くと、遠くの別の一羽がギーギーとこれまた耳障り悪く啼き返して、20分くらい交互にラップバトルを戦わせます。

春の早朝に布団の中で小鳥が啼きあっているのを聞くと、かわりばんこに音を奏でていて下記の曲のごとく心地よいものですが、冬の真昼間にラップバトルはちょっとほどほどに遠慮してもらいたいw


バッハ ブランデンブルク協奏曲第5番 第一楽章、アレグロ

【寝言】地震かな?

星の子チョビンのフクロウ

星の子チョビンのフクロウの動画がうpられていないというこの世のわびしさ。

f:id:fortran66:20210214005543j:plain

ふくろう 声 - 八奈見乗児 いつも木の枝に止まっているが、衝撃やチョビンの泣き声(物凄い強烈な音)などで枝から落ちて「地震かなァ」という(毎回そのシーンしか無い)。 星の子チョビン - Wikipedia

https://storage.tenki.jp/storage/static-images/forecaster_diary/image/1/11/114/11462/main/20210213232153/large.jpg

tenki.jp

f:id:fortran66:20210214012948p:plain

対数でプロットして欲しい。 f:id:fortran66:20210214012905p:plain

昨日は新月


星の子チョビン


藍美代子♪8.星のしずくの子守唄 Miyoko Ai オフィシャルYouTube チャンネル