元ネタはこのFortran プログラム。さして長いわけではないので全体を再現しておく。
c123456789012345678901234567890123456789012345678901234567890123456789012
c ----- -----
c ----- Vertex CTMRG for Ising model: Chess Board -----
c ----- version. -----
c ----- -----
c ----- '2002/08/24, Ver.1.0 by T.Nishino -----
c ----- -----
module var_size
integer,parameter :: mIT = 10000 ! Maximum Iteration Number
integer,parameter :: iq = 2 ! Ising (2-state) spins
integer,parameter :: mXX = 80 ! Maximum value of 'm'
integer,parameter :: mBW = iq**4 ! Size of the Boltzmann Weight
integer,parameter :: mXT = iq**3
& *mXX**2 ! Variable Size
end module
c
program main; use var_size
integer L, LL ! System Size L, Max. System Size LL
integer m, mX ! Block Spin State m, Max. Block Spin State mX
integer mN ! The value of m in the NeXT iteration.
integer tp ! Temporary Variable
c
real(8) K ! Inverse Temperature 1/T
real(8) mg ! Spontaneous Magnetization
real(8) u ! Bond Energy
real(8) Z ! Normalized Partition Function
real(8) eP, fP ! Normalization Const. for HRTM
real(8) eC, fC ! Normalization Const. for CTM
real(8) F(0:mIT) ! Free Energy
c
real(8) W(0:mBW) ! Boltzmann Weight
real(8) C(0:mXT) ! Corner Transfer Matrix
real(8) P(0:mXT) ! Half Row Transfer Matrix
real(8) G1(0:mXT) ! Work Area
real(8) G2(0:mXT) ! Work Area
real(8) G3(0:mXT) ! Work Area
real(8) G4(0:mXT) ! Work Area
c
call getprm( L, mX, K ) ! Set Physical Parameters
call initbw( K, W ) ! Initialize Boltzmann Weights
c ! Set CTM, HRTM, etc.
LL = 2; eC = 0.0d0; C(0)= 1.0d0
m = 1; eP = 0.0d0; P(0)= 1.0d0; P(1)= 0.0d0
goto 20
c
10 continue
call rjoint( m, m, iq*m, C, P, G1 )
call extP( iq, m, W, P, G2 )
call extC( iq, m, G1, G2, G3 )
call sjoint( iq *m, iq*m, G3, C )
call sjoint( iq *m, iq*m, C, G1 )
call observ( iq, m, mg, Z, G1 )
c
call diag( iq*m, G1, C, P ) ! C/P_temp
mN = min( iq*m, mX )
tp = iq*m*(iq*m-mN )
call renoP( iq,iq*m, mN, G2, P, G1(tp) )
call renoC( iq*m, mN, G3, C, G1(tp) ); m = mN
F(LL)= eP*8 + eC*4 + log(Z)
c
if(LL.gt.8) then; write(6,100) LL,mg,(F(LL)+F(LL-4)-2*F(LL-2))/8
else; write(6,100) LL,mg,0.0d0
end if
100 format('#',I4, F16.10, F16.10 )
c
20 continue
call symchP( iq, m, P, G1 ) ! Check the Symmetry of HRTM
call symchC( m, C, G1 ) ! Check the Symmetry of CTM
call norm( iq *m**2, P, fP ) ! Normalize HRTM
call norm( m**2, C, fC ) ! Normalize CTM
eC= eC + 2*eP + log(fC) ! Normalization Const for CTM
eP= eP + log(fP) ! Normalization Const for HRTM
LL = LL + 2; if( LL .LT. L ) goto 10
end program main
ccccc
subroutine errstp( ch ) ; character(8) ch
write(6,*) 'errstp <',ch,'>'; stop "Input_Error"; end
ccccc
subroutine getprm( L, mX, K ); use var_size
integer L, mX; real(8) K, T
write(6,*)'Max. Iteration # L'; read(5,*) L ; if( L.GT.mIT )
& call errstp( 'L > mIT')
write(6,*)'Maximum Dim. mX>= m'; read(5,*) mX; if( mX.GT.mXX )
& call errstp( 'mX > mXX')
write(6,*)'Parameter T (J=1) '; read(5,*) T; K = 1.0d0 / T
write(6,*)' mX = ', mX, ' T = ',T
end
ccccc
subroutine initbw( fK, W ); use var_size
integer i,j,k,l; real(8) fK, W(0:*) ! i---j
do i = 0, iq-1; do j = 0, iq-1 ! | W |
do k = 0, iq-1; do l = 0, iq-1 ! l---k
W(((i*iq+j)*iq+k)*iq+l)= exp( 4*fK*(i-0.5d0)*(j-0.5d0)
& + 4*fK*(j-0.5d0)*(k-0.5d0)
& + 4*fK*(k-0.5d0)*(l-0.5d0)
& + 4*fK*(l-0.5d0)*(i-0.5d0) )
end do; end do; end do; end do; end
cccccc
subroutine exchg2( m1, m2, A, B )
integer i, j, m1, m2; real(8) A(0:*), B(0:*)
do i = 0, m1-1
do j = 0, m2-1; B(j*m1+i) = A(i*m2+j)
end do; end do; end
cccccc
subroutine exchg3( m1, m2, m3, A, B )
integer i, j, m1, m2, m3; real(8) A(0:*), B(0:*)
do i = 0, m1-1
do j = 0, m2-1; B((j*m1+i)*m3:(j*m1+i)*m3 + m3-1)
& = A((i*m2+j)*m3:(i*m2+j)*m3 + m3-1)
end do; end do; end
cccccc
subroutine observ( q, m, mg, Z, G )
integer i,j, q, m; real(8) mg, Z, G(0:*)
c
mg = 0.0d0 ; Z = 0.0d0
do i = 0, q-1
do j = 0, m-1
Z = Z + G(((i*m+j)*q+i)*m+j)
mg = mg + G(((i*m+j)*q+i)*m+j)*(2*i-1)
end do; end do; mg = mg / Z
end
cccccc
subroutine sjoint( m1, m2, A, B )
integer i,j, m1, m2; real(8) A(0:*), B(0:*)
B(0:m2*m2-1)= 0.0d0
do i = 0, m1-1
do j = 0, m2-1; B(j*m2:j*m2 + m2-1) = ! *m1* m2
& B(j*m2:j*m2 + m2-1) + A(i*m2+j) ! *m1* m2
& * A(i*m2:i*m2 + m2-1) ! * * m2 m2
end do; end do; end
cccccc
subroutine ujoint( m1, m2, m3, A, B, C )
integer i,j, m1, m2, m3; real(8) A(0:*), B(0:*), C(0:*)
C(0:m2*m3-1)= 0.0d0
do i = 0, m1-1
do j = 0, m2-1; C(j*m3:j*m3 + m3-1) = ! *m1* m2
& C(j*m3:j*m3 + m3-1) + A(i*m2+j) ! *m1* m3
& * B(i*m3:i*m3 + m3-1) ! * * m2 m3
end do; end do; end
cccccc - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
subroutine rjoint( m1, m2, m3, A, B, C )
integer m1, m2, m3
real(8) A(m2,m1), B(m3,m2), C(m3,m1)
C = MATMUL(B,A)
end
cccccc - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
subroutine extP( q, m, W, P, G )
integer q, m; real(8) W(0,*), P(0:*), G(0:*)
c
call exchg2( m, q*m, P, G ) ! iaj -> [a]ji
call ujoint( q, m*m, q**3, G, W, P ) ! [a] bcd
call exchg2( m,m *q**3, P, G ) ! ibcdj <---- jibcd
call exchg3( m, q,m *q**2, G, P ) ! bicdj
c ! i b
G(0:m**2*q**3-1) = P(0:m**2*q**3-1) ! +-a-+-c
end ! j d
cccccc
subroutine extC( q, m, G1, G2, G3 )
integer q, m; real(8) G1(0:*), G2(0:*), G3(0:*)
c
call exchg3( m, q, m, G1, G3 ) ! iaj -> [ai]j
call ujoint( q*m, m, q**2*m, G3, G2, G1 ) ! [ai] bck
call exchg3( m, q, q *m, G1, G3 ) ! bjck <--- jbck
end
ccccc
subroutine diag( N, A, E, tp )
integer N, info; real(8) A(0:*), E(0:*), tp(0:*)
c
call DSYEV('V','U', N, A, N, E, tp, (N+2)*N, info )
c-SX- call dcsmaa( A, N, N, E, tp, info )
if( info .ne. 0 ) then
write(6,*) 'dsyev_info = ', info; call errstp(' -diag- ')
end if
end
cccccc
subroutine renoP( q, qm, mN, G, P, R )
integer q, qm, mN; real(8) G(0:*), P(0:*), R(0:*)
call rjoint( mN, qm, q*qm, R, G, P ) ! x[ai]
call exchg3( mN, q, qm, P, G ) ! [ai]bcj -> xbcj
call exchg2( q*mN, qm, G, P ) ! bxcj
call rjoint( mN, qm, q*mN, R, P, G ) ! y[cj] <- cjbx
c ! [cj]bx
P(0:q*mN**2-1) = G(0:q*mN**2-1) ! y bx
end
ccccc
subroutine renoC( qm, mN, G, C, R )
integer qm, mN; real(8) G(0:*), C(0:*), R(0:*)
call rjoint( mN, qm, qm, R, G, C ) ! x[ai] y[bj]
call exchg2( mN, qm, C, G ) ! [ai]bj -> xbj -> [bj]x
call rjoint( mN, qm, mN, R, G, C ) ! y x
end
cccccc
subroutine norm( N, V, w )
integer N; real(8) V(0:*), w
call maxel( N, V, w ); V(0:N-1)= V(0:N-1)/w
end
ccccc
subroutine maxel( N, V, w )
integer i, N; real(8) V(0:*), w; w = 0.0d0
do i = 0, N-1; if( abs(V(i)) .GT. w ) w = abs(V(i))
end do; end
ccccc
subroutine symchP( q, m, P, tp )
integer i,j,k, q, m; real(8) P(0:*), tp(0:*), T, S; S = 0.0d0
do i = 0, m-1
do j = 0, q-1
do k = 0, m-1 ; tp((i*q+j)*m+k)
& = ( P((i*q+j)*m+k) + P((k*q+j)*m+i) )/2
if( S .lt. ABS( P((i*q+j)*m+k) - P((k*q+j)*m+i) ) )
& S = ABS( P((i*q+j)*m+k) - P((k*q+j)*m+i) )
end do; end do; end do
c
P(0:m*q*m-1) = tp(0:m*q*m-1); call maxel( m*q*m, P, T)
c
if( S/T .gt. 1.0E-12 ) write(6,*) 'sym P = ', S/T, q, m
end
cccccc
subroutine symchC( m, C, tp )
integer i,j, m; real(8) C(0:*), tp(0:*), T, S; S = 0.0d0
do i = 0, m-1
do j = 0, m-1 ; tp(i*m+j) =( C(i*m+j) + C(j*m+i) )/2
if( S .lt. ABS( C(i*m+j) - C(j*m+i) ) )
& S = ABS( C(i*m+j) - C(j*m+i) )
end do; end do
c
call maxel( m*m, C, T); C(0:m*m-1)= tp(0:m*m-1)
c
if( S/T .gt. 1.0E-12 ) write(6,*) 'sym C = ', S/T, m
end
このプログラムを見てとりあえず疑問に思ってほしいことは、
多次元配列を使わないのは何故だ？
ということなわけだ。 Fortran の C に比べた(K&R C に比べた)利点の一つは、多次元配列が数値計 算で使いたいようなやりかたで使えるということである。いわゆる整合配列というものね。
もっとも、C99 ではこれは標準の機能としてサポートされるようになった。 複素数型がないことや総称関数がないといった点についても Fortran とほぼ同じ形でのサポートが追加されたので新しくプログラム開発を するのに Fortran を使う必要は、どうしても Fortran でないと性能がでないタ コな機械を使わざると得ないといった悲惨な場合を除いてはほぼ無くなったといっ ていいと思う。まあ、結構そういう場合はあるんだけどさ。
と、問題はそういうことではなく、何故多次元配列を使わないのかということ である。 本人の主張では、「自分で見易い」ためとのことなので、きっと本人には それが見易いんだろうしデバッグもきっと本人はそれが楽なんだと思われる。
が、私としては、これがわかりやすいプログラムであると思うようなセンスは 嫌いだ。以下、サブルーチン一つだけを例に何故嫌いかを議論していこう。
subroutine symchC( m, C, tp )
integer i,j, m; real(8) C(0:*), tp(0:*), T, S; S = 0.0d0
do i = 0, m-1
do j = 0, m-1 ; tp(i*m+j) =( C(i*m+j) + C(j*m+i) )/2
if( S .lt. ABS( C(i*m+j) - C(j*m+i) ) )
& S = ABS( C(i*m+j) - C(j*m+i) )
end do; end do
c
call maxel( m*m, C, T); C(0:m*m-1)= tp(0:m*m-1)
c
if( S/T .gt. 1.0E-12 ) write(6,*) 'sym C = ', S/T, m
end
が元のプログラムである。これをもう少し普通の形に直すと、以下のようにな
るだろう。
subroutine symchC(m,C,tp)
integer i,j, m
real*8 C(m,m),tp(m,m), T, S
call maxel( m*m, C, T)
S = 0.0d0
do i = 1, m
do j = 1, m
tp(i,j) =(C(i,j) + C(j,i) )*0.5D0
S =MAX(S,ABS(C(i,j)-C(j,i)))
end do
end do
c
C(1:m*m)= tp(1:m*m)
c
if( S/T .gt. 1.0E-12 ) write(6,*) 'sym C = ', S/T, m
end
maxel を呼ぶ場所が変わっているが、これは(多分)結果に影響しない。という
か、これは配列の要素の絶対値の最大値を求めるだけなので、 C を書換える
前ならいつ呼んでも同じと思われる。
書き直したもののほうなら、
まあ、インデントは流儀があるのは確かだが、 一つ重要な点として、「機械(エディタ)が自動的に付けられるものであるこ とが必須」である。これにより、人間がインデントを正しく付けるという苦労 から解放されるだけでなく、ネスティングのエラーをエディタだけで発見でき るというおまけもついてくる。
また、書き直したものでは、最大値を求めるのに if 文を使わないで max 組み込み関数を使っていることに注意して欲しい。こうすることで ABS(C(i,j)-C(j,i)) がプログラム中に2回発生するのを回避できる。
同じ表現が複数回でてくることは以下の理由で好ましくない
ここで重要なのは、もちろん最初のほう、つまり、人間の間違いの可能性を増 やすということである。性能的には、コンパイラによっては MAX を使うほう が遅いということもありえる。この場合には、中間変数を使ってでも表現の重 複を避けることが望ましいと私は思う。つまり
tmp = ABS(C(i,j)-C(j,i)) if (tmp .gt. s) s = tmpとするのである。元のプログラム(をとりあえず添字だけ直した)
if( S .lt. ABS(C(j,i)-C(i,j)))
& S = ABS(C(j,i)-C(i,j))
に比べると文字数でも結局減っている。一般にプログラムのバグの数は文字数
(というかトークン数)に比例すると言われており、より簡潔にプログラムを表
現することでよりバグが少なくまた発見もしやすいプログラムを書くことがで
きるのである。
さて、対称行列を作るなら、配列を2つ用意しなくても以下のようにすれば十
分である。
subroutine symchC(m,C)
integer i,j, m
real*8 C(m,m), T, S
call maxel( m*m, C, T)
S = 0.0d0
do i = 1, m-1
do j = i+1, m
S =MAX(S,ABS(C(i,j)-C(j,i)))
C(i,j) = (C(i,j) + C(j,i))*0.5D0
C(j,i) = C(i,j)
end do
end do
c
if( S/T .gt. 1.0E-12 ) write(6,*) 'sym C = ', S/T, m
end
つまり、配列の上(か下かは定義の問題だが)半三角だけをカバーするループを
回して、平均値を計算したら上と下の両方に代入する。
この形では、余計な配列を必要としない分このサブルーチンの独立性が高くなっ ており、その分バグの可能性が減っている。もっとも、2重ループの範囲を正 しく設定するのに多少の注意が必要なことも確かであり、 書く時にバグが入 りやすいという問題が新しく導入されてしまっている。
しかしながら、元の形では、「何故同じ計算を2回して必要もない配列を使っ て結果をしまうのか」という疑問も読む人によっては発生するであろう。これ に対して最終版ではそのような心配はとりあえずない。
もっとも、速度が問題ならば最終版のほうが明らかに望ましい。計算量が半分 以下であり、その結果計算時間が半分程度になるはずだからである。
というわけで私の主張は