fortran66のブログ

fortran について書きます。

【メモ帳】暇つぶしに Julia で Camplbell-Baker-Hausdorff 計算してみる

Camplbell-Baker-Hausdorff の公式

Camplbell-Baker-Hausdorff の公式という割に、実際の展開式を導いたのはソ連の Dynkin らしいので、この名でいいのかとも思いますが、アカと間違えられて FBI に通報されても困るので、長い物には巻かれましょうw

exp(A)exp(B) = exp(C) の時、C=A+B+[A, B]/2 + ([A,[A,B]]-[B,[A,B]])/12 + [B,[A,[A, B]]]/24+...

Julia 言語での行列積・行列 exp

Julia 言語の行列の積 X*Y は Fortran での matmul(X, Y) になっており、exp(X) も数学的な定義に従って exp のべき展開での行列積になっており、Fortran での行列要素の exp(x_ij) ではないので、これを利用してCamplbell-Baker-Hausdorff の公式に適当な行列を入れて計算して見ることにします。

A=[1 1 
   1 1]
2×2 Array{Int64,2}:
 1  1
 1  1
B = [2 1
     1 2]/10
2×2 Array{Float64,2}:
 0.2  0.1
 0.1  0.2
f(x, y)=x*y-y*x
f (generic function with 1 method)
f(A, B)
2×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0
exp(A)* exp(B)
2×2 Array{Float64,2}:
 5.53968  4.43451
 4.43451  5.53968
exp(A+B+f(A,B)/2)
2×2 Array{Float64,2}:
 5.53968  4.43451
 4.43451  5.53968
A=[1 1 
   1 0]/10
2×2 Array{Float64,2}:
 0.1  0.1
 0.1  0.0
f(A, B)
2×2 Array{Float64,2}:
  0.0   0.01
 -0.01  0.0 
exp(A)* exp(B)
2×2 Array{Float64,2}:
 1.37607   0.26518
 0.252291  1.24676
exp(A+B+f(A,B)/2+(f(A, f(A, B))-f(B, f(A, B)))/12)
2×2 Array{Float64,2}:
 1.37607  0.265202
 0.25227  1.24676 
exp(A+B+f(A,B)/2+(f(A, f(A, B))-f(B, f(A, B)))/12 -(f(B, f(A, f(A, B))))/24  )
2×2 Array{Float64,2}:
 1.37607   0.26518
 0.252292  1.24676

行列 A, B が可換の時は、当然ながら C=A+B で exp(A)exp(B)=exp(A+B) が成り立ちます。

非可換の時は A, B を 0.1 のオーダーにしたので、交換子は 0.01 のオーダーで、四次の項まで取ることで小数点以下5桁くらいまで一致するようです。まぁそんなものかなという気のする結果ではあります。

Julia 言語では 行列積の関数も適当に書けるので、確かに楽だなと思います。プロトタイプ向きなのがよく分かります。

【メモ帳】Fortran 2018 C言語の相互運用 MFE §21.6 その2

MFE §21.6.5 Allocatable objects

C 側で Fortran の allocatable 変数を deallocate/allocate します。

これを使えば、ZMQ の zmq_msg_init_data も出来ると思います。

fortran66.hatenablog.com

実行結果

gfortran-9/gcc-9

hp8@HP8:~/f2018$ gfortran-9 mfe2165.c mfe2165.f90
hp8@HP8:~/f2018$ ./a.out
 GCC version 9.1.0
 -mtune=generic -march=x86-64

 x
allocatable
rank=2
lowerbound1 extent5 spacing4
allocated
dealloca retcode 0
unallocated
alloca retcode 0
 T           4           5
           0           1
           3           5
   1.00000000       1.00000000       1.00000000       1.00000000       1.00000000       1.00000000       1.00000000       1.00000000       1.00000000       1.00000000       1.00000000       1.00000000       1.00000000       1.00000000       1.00000000       1.00000000       1.00000000       1.00000000       1.00000000       1.00000000

icc/ifort

hp8@HP8:~/f2018$ icc -c mfe2165.c
hp8@HP8:~/f2018$ ifort mfe2165.o mfe2165.f90
hp8@HP8:~/f2018$ ./a.out
 Intel(R) Fortran Intel(R) 64 Compiler for applications running on Intel(R) 64,
 Version 19.1.0.056 Pre-Release Beta Build 20190321


 x
allocatable
rank=2
lowerbound1 extent5 spacing4
allocated
dealloca retcode 0
unallocated
alloca retcode 0
 T           4           5
           0           1
           3           5
   1.000000       1.000000       1.000000       1.000000       1.000000
   1.000000       1.000000       1.000000       1.000000       1.000000
   1.000000       1.000000       1.000000       1.000000       1.000000
   1.000000       1.000000       1.000000       1.000000       1.000000

ソース・プログラム

program test
    use, intrinsic :: iso_fortran_env
    use, intrinsic :: iso_c_binding
    implicit none
    interface
        subroutine alloc(cdesc) bind(c, name = 'alloc')
            use, intrinsic :: iso_c_binding
            real, intent(in out), allocatable :: cdesc(..)
        end subroutine alloc
    end interface

    real, allocatable :: x(:, :)

    print *, compiler_version()
    print *, compiler_options()
    print *

    print *, 'x'

    allocate(x(5, 5))
    call alloc(x)
    print *, allocated(x), shape(x)
    print *, lbound(x)
    print *, ubound(x)
    x = 1.0
    print *, x
    deallocate(x)

end program test
#include <stdio.h>
#include "ISO_Fortran_binding.h"

int alloc(CFI_cdesc_t *cdesc) {
    CFI_attribute_t attrib = cdesc->attribute;
    if (attrib == CFI_attribute_allocatable) {printf("allocatable\n");};
    if (attrib == CFI_attribute_pointer    ) {printf("pointer\n");};
    if (attrib == CFI_attribute_other      ) {printf("other\n");};



    CFI_rank_t nrank = cdesc->rank;
    printf("rank=%d\n", (int)nrank);

    CFI_dim_t *dim = cdesc->dim;
    CFI_index_t lowerbound = dim->lower_bound;
    CFI_index_t extent     = dim->extent;
    CFI_index_t sm         = dim->sm;
    printf("lowerbound%d extent%d spacing%d\n", (int) lowerbound, (int) extent, (int) sm);


    if (cdesc->base_addr) {
        fputs("allocated\n", stderr);
        int iret = CFI_deallocate(cdesc);
        printf("dealloca retcode %d\n", iret);
    }

    fputs("unallocated\n", stderr);
    CFI_index_t lo[] = {0, 1};
    CFI_index_t hi[] = {3, 5};
    size_t len = 0;
    int iret = CFI_allocate(cdesc, lo, hi, len);
    printf("alloca retcode %d\n", iret);

}

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)

【メモ帳】○○的の「的」について

杉山孚富 著「漢文故事熟語詳解」より

明治36年 (1903年)初版 明治38年(1905年)訂正6版

面倒なので当用漢字でw カタカナもひらがなに。その他適当に濁点づけ。

てノ部 的 「近来人好んで的の字を用ゆ何々的など云へり此は元来助字にして意味あるにあらず句調を助ける添文字なり晋人好んで馨の字用ひ宋人好んで底の字を用ひ清人好んで兒の字を用ゆるに同じ皆な助字なり音学五書に云う的の字入声となれば則ち「チャク」となり去声に転ずれば「テフ」となり上声に転ずれば「テイ」となる宋人の書中に凡ての語助は皆底(テイ)に作る的の字なしと然らば則ち的も底も同字として助字に用ゆること知るべし」

角川「新字源」には英語の ...tic の意とありますが、19世紀末から20世紀初頭に流行り始めた用法なのでしょうか? 近時支那では助詞の「の」とか of みたいな感じな気がしますが。

「漢文故事熟語詳解」

漢文故事熟語詳解 (金刺芳流堂): 1902|書誌詳細|国立国会図書館サーチ

国会図書館データの年号がずれていて謎。 奥付の方には「受験応用 漢文故事熟語詳解」と。学研の付録みたいな雑然としたオモシロ本。

【寝言】ペンス副大統領の反中演説ようやく!?

ペンス演説第二段

当初、天安門事件30周年の日、六月四日にする予定で、その後延び延びになっていた、マイク・ペンス米副大統領の反中演説が、十月一日の中共建国記念日に向けて行われるようです。

www.youtube.com

去年十月にもハドソン研究所で画期的な演説をしていますので、期待が持てます。はやく完全なる封じ込めを宣言して欲しいです。

fortran66.hatenablog.com

fortran66.hatenablog.com

www.sankei.com

f:id:fortran66:20190907232343j:plain
マイク・ペンス

 【ワシントン=黒瀬悦成】米ホワイトハウスのミラー副大統領副報道官は6日、ペンス副大統領が中国に関する政策演説を今秋、ワシントン市内の政策研究機関「ウィルソン・センター」で行うことを明らかにした。中国の習近平体制による新疆ウイグル自治区ででの人権抑圧や香港情勢などで習体制を批判する内容になるとみられ、貿易や安全保障に加え、人権・民主化問題でも中国に全面的圧力を加えていく立場を鮮明に打ち出す考えだ。

 ペンス氏の演説は当初、6月に予定されていたが、同月末にトランプ大統領と中国の習近平国家主席との首脳会談を控え、米中の貿易交渉が一つの山場を迎えていたため、交渉進展を優先させる思惑から延期された経緯がある。

 ペンス氏は昨年10月にも政策研究機関「ハドソン研究所」で、米政権による中国との「全面対決」を宣言する演説を行っており、今度の演説はその「第2弾」に位置づけられる。

 ペンス氏は、最近では香港情勢に関し、中国が米国と貿易問題で合意したいのであれば、1984年の中英共同宣言に基づく香港の自治権を尊重すべきだと指摘し、中国が抗議デモを武力鎮圧すれば「合意は困難となる」と警告した。

 一方、ポンペオ国務長官は6日、中西部カンザス州の大学で講演し、今月中旬からニューヨークで始まる国連総会の場で、中国による新疆ウイグル自治区での住民弾圧に関し、各国に「中国糾弾」を呼びかけていく考えを明らかにした。

 ポンペオ氏は、中国当局自治区住民らに対する処遇は「今世紀の世界における最悪の汚点になる可能性がある」と指摘し、「米国は彼らに自由がもたらされることを望む。ことは彼らの基本的かつ誰にも奪うことのできない権利に関わる問題だ」と強調した。

【メモ帳】Fortran 2018 C言語の相互運用 MFE §21.6

MFE §21.6 Accessing Fortran objects

~21.6.4

実行結果

gfortran-9

hp8@HP8:~/f2018$ gfortran-9 c_mfe216.c mfe216.f90
hp8@HP8:~/f2018$ ./a.out
 GCC version 9.1.0
 -mtune=generic -march=x86-64

 x
contiguous
associated

 ip => null()
discontiguous
disassociated

 ip => m(::2)
discontiguous
associated

ifort

hp8@HP8:~/f2018$ icc -c c_mfe216.c
hp8@HP8:~/f2018$ ifort c_mfe216.o mfe216.f90
hp8@HP8:~/f2018$ ./a.out
 Intel(R) Fortran Intel(R) 64 Compiler for applications running on Intel(R) 64,
 Version 19.1.0.056 Pre-Release Beta Build 20190321
 -traceback

 x
contiguous
associated

 ip => null()
discontiguous
disassociated

 ip => m(::2)
discontiguous
associated

ソース・プログラム

Fortran

program test
    use, intrinsic :: iso_fortran_env
    use, intrinsic :: iso_c_binding
    implicit none
    interface
        subroutine iscont(cdesc) bind(c, name = 'iscont')
            use, intrinsic :: iso_c_binding
            type(*), intent(in) :: cdesc(..)
        end subroutine iscont
    end interface

    real :: x(10) = 1.0
    integer, target  :: m(10) = 8
    integer, pointer :: ip(:) => null()

    print *, compiler_version()
    print *, compiler_options()
    print *

    print *, 'x'
    call iscont(x)

    print *, 'ip => null()'
    call iscont(ip)

    print *, 'ip => m(::2)'
    ip => m(::2)
    call iscont(ip)

end program test

C

#include <stdio.h>
#include "ISO_Fortran_binding.h"

int iscont(CFI_cdesc_t *cdesc) {
    if (CFI_is_contiguous(cdesc)) {
        fputs("contiguous\n", stderr);
    } else {
        fputs("discontiguous\n", stderr);
    }

    if (cdesc->base_addr) {
        fputs("associated\n\n", stderr);
    } else {
        fputs("disassociated\n\n", stderr);
    }
}

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)

【メモ帳】Fortran 2018 C言語の相互運用 MFE §21.5

Modern Fortran Explained - 2018 §21.5 C descriptors

暑いので適当w interface に type(*) :: x(..) を使ったのが新奇。 

attribute がうまく取れていないのが謎。

実行結果

gfortran-9/gcc-9

hp8@HP8:~/f2018$ gfortran-9 c_mfe215b.c mfe215b.f90
hp8@HP8:~/f2018$ ./a.out
 -mtune=generic -march=x86-64

other
rank=0
lowerbound0 extent0 spacing0
CFI_type_double

other
rank=2
lowerbound0 extent10 spacing4
CFI_type_float

other
rank=0
lowerbound0 extent10 spacing4
CFI_type_double

other
rank=1
lowerbound0 extent100 spacing4
CFI_type_float

other
rank=0
lowerbound0 extent100 spacing4
CFI_type_float

STOP normal termination

ifort/icc

hp8@HP8:~/f2018$ vi mfe215b.f90
hp8@HP8:~/f2018$ icc -c   c_mfe215b.c
hp8@HP8:~/f2018$ ifort   c_mfe215b.o mfe215b.f90
hp8@HP8:~/f2018$ ./a.out


other
rank=0
lowerbound13 extent1886220131 spacing1
CFI_type_double

other
rank=2
lowerbound0 extent10 spacing4
CFI_type_float

other
rank=0
lowerbound469778721 extent-604196248 spacing26577600
CFI_type_double

other
rank=1
lowerbound0 extent100 spacing4
CFI_type_float

other
rank=0
lowerbound0 extent0 spacing0
CFI_type_float

normal termination

ソース・プログラム

fortran

program test
    use, intrinsic :: iso_fortran_env
    use, intrinsic :: iso_c_binding
    implicit none
    interface
        subroutine cfi_members(string) bind(c, name = 'cfi_members')
            use, intrinsic :: iso_c_binding
            type(*) :: string(..)
        end subroutine cfi_members
    end interface

    real(kind=kind(1.0d0)), target :: a
    real :: b(10, 10)
    real(kind=kind(1.0d0)), pointer :: p
    real, allocatable :: x(:), y
    integer :: i, j

    print *, compiler_options()
    print *

    a = 1.0
    call cfi_members(a)

    forall(i=1:10, j=1:10) b(i, j) = i + j
    call cfi_members(b)

    p => a
    call cfi_members(p)

    x = [(i, i = 1, 100)]
    call cfi_members(x)

    call cfi_members(y)

    stop 'normal termination'
end program test

C

#include <stdio.h>
#include "ISO_Fortran_binding.h"

void cfi_members(CFI_cdesc_t *string) {
    CFI_attribute_t attrib = string->attribute;
    if (attrib == CFI_attribute_allocatable) {printf("allocatable\n");};
    if (attrib == CFI_attribute_pointer    ) {printf("pointer\n");};
    if (attrib == CFI_attribute_other      ) {printf("other\n");};

    CFI_rank_t nrank = string->rank;
    printf("rank=%d\n", (int)nrank);

    CFI_dim_t *dim = string->dim;
    CFI_index_t lowerbound = dim->lower_bound;
    CFI_index_t extent     = dim->extent;
    CFI_index_t sm         = dim->sm;
    printf("lowerbound%d extent%d spacing%d\n", (int) lowerbound, (int) extent, (int) sm);

    CFI_type_t ctype = string->type;
    switch (ctype) {
        case CFI_type_int:
//        case CFI_type_int32_t:
            printf("CFI_type_int\n");
            break;
        case CFI_type_Bool:
            printf("CFI_type_Bool\n");
            break;
        case CFI_type_float:
            printf("CFI_type_float\n");
            break;
        case CFI_type_double:
//        case CFI_type_long_double:
            printf("CFI_type_double\n");
            break;
        case CFI_type_cptr:
            printf("CFI_type_cptr\n");
            break;
//        case CFI_type_int64_t:
//        case CFI_type_long:
//        case CFI_type_long_long:
//        case CFI_type_size_t:
//        case CFI_type_intmax_t:
//        case CFI_type_intptr_t:
//          case CFI_type_ptrdiff_t:

        case CFI_type_short:
//        case CFI_type_int16_t:
            printf("CFI_type_short\n");
            break;
        case CFI_type_int8_t:
//        case CFI_type_signed_char:
            printf("CFI_type_signed_char\n");
            break;
// undefined        case CFI_type_least8_t:
//        case CFI_type_least16_t:
//        case CFI_type_least32_t:
//        case CFI_type_least64_t:
//        case CFI_type_fast8_t:
//        case CFI_type_fast16_t:
//        case CFI_type_fast32_t:
//        case CFI_type_fast64_t:
        case CFI_type_struct:
            printf("CFI_type_struct\n");
            break;
        case CFI_type_other:
            printf("CFI_type_other\n");
            break;
        case CFI_type_float_Complex:
            printf("CFI_type_float_Complex\n");
            break;
        case CFI_type_double_Complex:
//        case CFI_type_long_double_Complex:
            printf("CFI_type_double_Complex\n");
            break;
        case CFI_type_char:
            printf("CFI_type_char\n");
            break;
        default:
            printf("case default\n");
    }



    printf("\n");
}

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)

【寝言】覗き魔 google 通販履歴貯めてるw その他

Googleの「自分のデータをダウンロード」

お盆休みに、google の「自分のデータをダウンロード」を、ほんの一部だけ見てみたのですが、通販での買い物履歴が json データになってたまっていてムカつきました。amazon などは別メアドなのに、スパムフィルター代わりに gmail を通していたら、メールの中身を読んでいて amazon の買い物記録もバッチシまとめられていました。ヤマジュンのホモマンガとか。

labaq.com

ウホッ!!いい男たち~ヤマジュン・パーフェクト

ウホッ!!いい男たち~ヤマジュン・パーフェクト

緊急地震速報

緊急地震速報の「ギュッ!ギュッ! 地震が来ます。 ギュッ!ギュッ! 地震が来ます。」というのを聞くと、いつも魔法陣グルグルの魔物を知らせる人形を思い出します。

f:id:fortran66:20190907012959p:plain
魔法陣グルグル

EMOTION the Best 魔法陣グルグル DVD-BOX 1

EMOTION the Best 魔法陣グルグル DVD-BOX 1

月と木星

今日は上弦くらいの月の下のすぐ近くに木星が見えていました。とてもきれいでした。