IDO Memo - @zigen 's note

非保存型IDO法の説明

IDO-NCF
the no conservative form of the Interpolated Differential Operator (IDO- CF)

今日(2007.3.1),東京工業大学を訪問し、スパコンツバメなどを見学。
ツバメは一般学生もIDを作れて特別な許可なしで32CPUまで使えるとのこと…すげー
また、SUNの講習会があったとかでいろいろ話を聞けました

1. ツバメは各CPU毎にメモリが割り与えられているためそれがボトルネックになったりする
2. キャッシュをうまく使うことによって８time（8倍）の性能差が現れたりする
3. f(i,j)等をループする場合下記の違いでさえかなりの性能差が出る

      Do 100 i=min, max
      Do 100 j=min, max
100   continue

or

      Do 100 j=min, max
      Do 100 i=min, max
100   continue

こんなの全部コンパイラがチューニングしてくれるのが一番いいんだろうけど、やっぱ奥が深いｗ
Solaris使ってみてはどうか？とのこと、Solarisには前々から興味があったのだが、コンパイラもSUN純正のがついてくるらしいし、iccとかも学生はただらしい（ぜんぜん知りませんでした）。それにクラスターのジョブ管理は管理用のアプリ（ウェブブラウザ上で動く）のがあるらしくて、各ノードがどれぐらい稼動してるか簡単に見れたのがすごかった。

1. スイッチについてもGigaEntherじゃなくてもっといいのがあると言ってたな（今度調べとけ）

ここからが本題Vlasov方程式を解くならIDO使ってみてはとのこと、（位相空間での）1次元（xとvでの二次元）ではIDOで解けない条件はないだろうとまで言われてた、今度IDOについても調べとけ

テストコード

ido_rk4.f90

     Program main
     implicit none
     double precision :: dt, dx, o_dx, o_dx2, v_v, f_d, dx_isign, xx
     double precision :: k1, k2, k3, k4, f1, f2, f3, f4
     double precision :: kx1, kx2, kx3, kx4, fx1, fx2, fx3, fx4
     double precision :: dft, dftx
     double precision, dimension(:), allocatable :: f, fn, g, gn, ft, ftx, initial
     double precision, dimension(:), allocatable :: f1_tmp, f1x_tmp, f2_tmp, f2x_tmp
     double precision, dimension(:), allocatable :: f1t_tmp, f1tx_tmp, f2t_tmp, f2tx_tmp
     double precision, dimension(:), allocatable :: f3t_tmp, f3tx_tmp, f4t_tmp, f4tx_tmp
  integer :: i, j, t, x_up, n_x = 20

     allocate( initial(0:n_x+1), f(0:n_x+1), fn(0:n_x+1), g(0:n_x+1), gn(0:n_x+1), ft(0:n_x+1), ftx(0:n_x+1) )
     allocate( f1_tmp(0:n_x+1), f1x_tmp(0:n_x+1), f2_tmp(0:n_x+1), f2x_tmp(0:n_x+1) )
     allocate( f1t_tmp(0:n_x+1), f1tx_tmp(0:n_x+1), f2t_tmp(0:n_x+1), f2tx_tmp(0:n_x+1) )
     allocate( f3t_tmp(0:n_x+1), f3tx_tmp(0:n_x+1), f4t_tmp(0:n_x+1), f4tx_tmp(0:n_x+1) )

!=======================================================================
! Initial Setting
!=======================================================================

     dx = 0.01d0
     o_dx = 1.d0/dx
     dt = 0.01d0
  dx_isign = 1.d0
     v_v = 0.1d0
     xx = -v_v*dt
     dx_isign = -dsign(1.d0, v_v)

     do i = 0, n_x+1
       initial(i) = 0.d0
       f(i) = 0.d0
	fn(i) = 0.d0
       ft(i) = 0.d0
	ftx(i) = 0.d0
	f1_tmp(i) = 0.d0
	f1x_tmp(i) = 0.d0
	f2_tmp(i) = 0.d0
	f2x_tmp(i) = 0.d0
       gn(i) = 0.d0
	  if( 5<=i .and. i<=10 ) f(i) = 1.d0
	  if( 5<=i .and. i<=10 ) initial(i) = 1.d0
     enddo
     do i = 0, n_x
       g(i) = (f(i) + f(i+1))/dx
     enddo
       g(n_x+1) = 0.d0

!=======================================================================
! Main Loop
!=======================================================================

     do t = 1, 10

!########################################################################

       do i = 1, n_x
	  x_up = i + int(dx_isign)
	  call cip(f(i), f(x_up), g(i), g(x_up), v_v, ft(i), ftx(i))
	  f1t_tmp(i) = ft(i)
	  f1tx_tmp(i) = ftx(i)

! f1_tmp(i) = dt*ft(i) !k1 = h*f(x_n,y_n)
! f1x_tmp(i) = dt*ftx(i)

       enddo

!########################################################################

       do i = 1, n_x
	  x_up = i + int(dx_isign)
	  call cip(f(i)+f1t_tmp(i)*dt/2.d0, f(x_up)+f1t_tmp(x_up)*dt/2.d0,  &
                  g(i)+f1tx_tmp(i)*dt/2.d0, g(x_up)+f1tx_tmp(x_up)*dt/2.d0,  &
                  v_v, ft(i), ftx(i))
         f2t_tmp(i) = ft(i) !k2 = h*f(x_n+dx,y_n+k1*dt)
         f2tx_tmp(i) = ftx(i)
       enddo

!########################################################################

       do i = 1, n_x
	  x_up = i + int(dx_isign)
	  call cip(f(i)+f2t_tmp(i)*dt/2.d0, f(x_up)+f2t_tmp(x_up)*dt/2.d0,  &
                  g(i)+f2tx_tmp(i)*dt/2.d0, g(x_up)+f2tx_tmp(x_up)*dt/2.d0,  &
                  v_v, ft(i), ftx(i))
         f3t_tmp(i) = ft(i) !k3 = h*f(x_n+dx,y_n+k1*dt)
         f3tx_tmp(i) = ftx(i)
       enddo

!########################################################################

       do i = 1, n_x
	  x_up = i + int(dx_isign)
	  call cip(f(i)+f3t_tmp(i)*dt, f(x_up)+f3t_tmp(x_up)*dt,  &
                  g(i)+f3tx_tmp(i)*dt, g(x_up)+f3tx_tmp(x_up)*dt,  &
                  v_v, ft(i), ftx(i))
         f4t_tmp(i) = ft(i)
         f4tx_tmp(i) = ftx(i)
       enddo

!########################################################################

	do i = 1, n_x
         f(i) = f(i) + (f1t_tmp(i) + 2.d0*f2t_tmp(i) + 2.d0*f3t_tmp(i) + f4t_tmp(i))*dt/6.d0
         g(i) = g(i) + (f1tx_tmp(i) + 2.d0*f2tx_tmp(i) + 2.d0*f3tx_tmp(i) + f4tx_tmp(i))*dt/6.d0
	enddo
     enddo

     open(99, file="result_ido_rk4.txt")
     do i = 1, n_x
       write(99,*) i, initial(i), f(i)
       write(*,*) i, initial(i), f(i)
     enddo
     close(99)
     End program main

!-----------------------------------------------------------------------

     subroutine cip(f_i, f_up, g_i, g_up, v_v, ft, ftx)
  implicit none
  double precision :: f_i, f_up, g_i, g_up, v_v, ft, ftx, dx, o_dx, dx_isign, b, f_d
       dx = 0.01d0
       o_dx = 1.d0/dx
       dx_isign = -dsign(1.d0, v_v)
	f_d = (f_up - f_i)*o_dx*dx_isign
	b = (3.d0*f_d - 2.d0*g_i - g_up)*o_dx*dx_isign
	ft = -v_v*g_i
	ftx = -2.d0*v_v*b
     return
  end subroutine cip

!-----------------------------------------------------------------------