秋分 23 日

Equinox を迎える喜び。国民の休日になっているのももっとも。χ を崇める。

DevCloud

intel の無料試用 cloud server で fortran が利用できますが、python 用の kernel しか用意されていません。

昔、MS Azure 用に試みた手法で fortran 用の kernel が動いたのでメモっておきます。

縦軸は絶対値を取った |zeta| の誤りですｗなお Pillow や sympy 等はインストールされていないので、コンソールから pip install *** する必要があります。

    program main
        use, intrinsic :: iso_fortran_env
        implicit none
        integer, parameter :: kd = real128
        integer :: i
        do i = 0, 1000
            write(10, *) i / 100.0, abs(zeta(cmplx(0.5_kd, i / 10.0_kd, kd)))
        end do

    contains

        complex(kd) function zeta(s) result(res)
            import, only : kd
            complex(kd), intent(in) :: s
            integer, parameter :: mmax = 100
            integer :: m, j
            complex(kd) :: c, d, e, w(mmax)
            real(kd), allocatable :: bino(:)
            bino = [1.0_kd]
            c = 1.0_kd / (1 - 2.0_kd**(1-s))
            d = (0.0_kd, 0.0_kd)
            do m = 1, mmax
                e = (0.0_kd, 0.0_kd)
                do j = 1, m
                    e = e - (-1)**j * bino(j) * j**(-s)
                end do
                d = d + e / 2.0_kd**m
                bino = [bino, 0.0_kd] + [0.0_kd, bino] ! binomial coefficients by Pascal's triangle
            end do
            res = c * d
        end function zeta

    end program main

%fig:
x, y1 = np.loadtxt('fort.10', unpack = True)
plt.grid(which='both')
plt.plot(x, y1)
plt.title('zeta function on 0.5 + x i', fontsize = 25)
plt.ylabel('|zeta|', fontsize = 20)

fortran66.hatenablog.com

coarray が簡単に試せるのがいいかと思います。gfortran しかない場合 OpenCoarray の準備が少しめんどくさい。

%fcflags: -coarray
    program Console3
        implicit none
        integer, parameter :: n = 10**6
        real :: x(n), y(n)
        
        integer :: me, np, ipi
        me = this_image()
        np = num_images()
        
        call random_init(.false., .true.) 
        call random_number(x)
        call random_number(y)
        ipi = count(x**2 + y**2 < 1.0)
        
        call co_sum(ipi, 1)
        if (me == 1) then 
            print *, 'number of images =', np, ' number of trials for Monte Carlo', n * np 
            print *, 4.0 * ipi / real(n * np) , 'error estimaton', 1.0 / (n * np)**0.5
        end if    
    end program Console3

fortran66.hatenablog.com

gfortran 用 kernel

github.com

from Terminal

git clone https://github.com/f66blog/azure

jupyter kernelspec install ~/azure/jupyter-gfort-kernel/gfort_spec --user

from Python3 cell

import sys

!{sys.executable} -m pip install -e ~/azure/jupyter-gfort-kernel --user

from terminal : edit kernel.json

which python

/blob/development-tools/versions/oneapi/2022.2/oneapi/intelpython/latest/bin/python

vi .local/share/jupyter/kernels/gfort_spec/kernel.json

"argv": [

" /blob/development-tools/versions/oneapi/2022.2/oneapi/intelpython/latest/bin/python",

intel fortran 用 kernel

github.com

from Terminal

git clone https://github.com/f66blog/jupyter-ifort-kernel.git

jupyter kernelspec install ~/jupyter-ifort-kernel/ifort_spec --user

from Python3 cell

import sys

!{sys.executable} -m pip install -e ../azure/jupyter-ifort-kernel --user

from terminal : edit kernel.json

which python

/blob/development-tools/versions/oneapi/2022.2/oneapi/intelpython/latest/bin/python

vi .local/share/jupyter/kernels/gfort_spec/kernel.json

"argv": [

" /blob/development-tools/versions/oneapi/2022.2/oneapi/intelpython/latest/bin/python",