あらゆるものを 1 次元で

2004/1/4 牧野


これは私の日誌の 2002/8 月辺りから。

元ネタはこの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 を書換える 前ならいつ呼んでも同じと思われる。 書き直したもののほうなら、 というのはまあ見ればわかると思う。元のプログラムでそれがわからないわけ ではないが、C(j,i)の代わりに C(i*m+j)と書くことでプ ログラムがわかりやすくなるとはとうてい思えない。また、インデント等に明 確な規則性が見えないこともプログラムの可読性を低めるのに貢献していると いえる。

まあ、インデントは流儀があるのは確かだが、 一つ重要な点として、「機械(エディタ)が自動的に付けられるものであるこ とが必須」である。これにより、人間がインデントを正しく付けるという苦労 から解放されるだけでなく、ネスティングのエラーをエディタだけで発見でき るというおまけもついてくる。

また、書き直したものでは、最大値を求めるのに 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回して必要もない配列を使っ て結果をしまうのか」という疑問も読む人によっては発生するであろう。これ に対して最終版ではそのような心配はとりあえずない。

もっとも、速度が問題ならば最終版のほうが明らかに望ましい。計算量が半分 以下であり、その結果計算時間が半分程度になるはずだからである。

というわけで私の主張は

というようなこと。