元ネタはこの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 endmaxel を呼ぶ場所が変わっているが、これは(多分)結果に影響しない。という か、これは配列の要素の絶対値の最大値を求めるだけなので、 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回して必要もない配列を使っ て結果をしまうのか」という疑問も読む人によっては発生するであろう。これ に対して最終版ではそのような心配はとりあえずない。
もっとも、速度が問題ならば最終版のほうが明らかに望ましい。計算量が半分 以下であり、その結果計算時間が半分程度になるはずだからである。
というわけで私の主張は