Gauss の掃き出し法で線型方程式を解く

#掃き出し法
よく教科書にあるようなものを使って、複素行列を複素型の変数で解く場合と、複素数と２ｘ２行列の同型構造を利用して、実数型の変数を使って二倍サイズの問題を解く場合を比較してみようと思い立ちました。

予想では、複素演算は実数演算の６倍くらいの演算量で、２倍サイズの線型方程式は計算量が n^3 に比例するので８倍くらいの演算量になって、これらに比例した計算時間になるだろうと思いました。

ところが、結果は２倍サイズ実数行列が劇遅で、複素行列が劇速となって、当惑ｗ

##寄り道
とりあえず、部分枢軸交換付きの Gauss 掃き出し法を使うので、素直に式通りにコーディングすると枢軸交換用の配列で配列の間接参照が出てしまうので、その効果を見てみます。

2000 x 2000 の乱数で生成した実行列と実ベクトルの Ax=b 型の線型方程式を解きます。コンパイラは Intel Fortran v.18 です。


###実行結果
結果を言うと、コンパイラが間接参照をうまく処理していて、一時変数に入れた場合と実行時間に差はありませんでした。以前は一時変数に入れてスカラー変数を配列の引数にしないと最適化が効きにくかったのですが、今はそういうことは無いようです。

```
time =    0.00    0.00:before gauss
time =   13.71   13.71:after  gauss
 RMS of diff.   2.9595977E-02
time =    0.01   13.72:before gauss
time =   13.48   27.20:after  gauss
 RMS of diff.   3.0429136E-02
```
ここでは並列化をかけて４コア８スレッドで３～４倍の高速化を成しています。

###プログラム
枢軸の交換に、ベクトルサブスクリプトを用いてオサレ感を出してみました。 

```fortran
ipiv([i, k]) = ipiv([k, i]) ! swap pivot
```
配列要素の交換にしか使えませんがｗ

```fortran
    module m_timing
      implicit none
    contains
      subroutine timimg(text)
        character(*), intent(in) :: text
        logical, save :: first = .true.
        integer, save :: ic0, ic1, ic2, icr
        if (first) then
          call system_clock(ic0, icr)
          ic1 = ic0
          first = .false.
        end if
        call system_clock(ic2, icr)
        print '(a, 2f8.2, 2a)', 'time =', real(ic2 - ic1) / icr, real(ic2 - ic0) / icr, ':', trim(text)
        ic1 = ic2
      end subroutine timimg
    end module m_timing
    
    module m_gauss
      implicit none
    contains
      subroutine gauss(a, b, x)
        real, intent(in out) :: a(:, :), b(:)
        real, intent(   out) :: x(:)
        integer :: i, j, k, ipiv(size(b))
        forall(i = 1:size(b)) ipiv(i) = i
        do i = 1, size(a, 2)
          k = maxloc(abs(a(ipiv(i:), i)), 1) + i - 1
          ipiv([i, k]) = ipiv([k, i]) ! swap pivot
          do j = i + 1, size(a, 1)     
            b(ipiv(j))    = b(ipiv(j))    - a(ipiv(j), i) / a(ipiv(i), i) * b(ipiv(i))
            a(ipiv(j), :) = a(ipiv(j), :) - a(ipiv(j), i) / a(ipiv(i), i) * a(ipiv(i), :) 
          end do
        end do
        do i = size(b), 1, -1
          x(i) = ( b(ipiv(i)) - sum(a(ipiv(i), i+1:) * x(i+1:)) ) / a(ipiv(i), i)
        end do  
      end subroutine gauss
      
      subroutine gauss0(a, b, x)
        real, intent(in out) :: a(:, :), b(:)
        real, intent(   out) :: x(:)
        real :: f
        integer :: i, j, k, m, ipiv(size(b))
        forall(i = 1:size(b)) ipiv(i) = i
        do i = 1, size(a, 2)
          k = maxloc(abs(a(ipiv(i:), i)), 1) + i - 1  
          ipiv([k, i]) = ipiv([i, k])
          m = ipiv(i)
          do j = i + 1, size(a, 1)
            k = ipiv(j)  
            f = a(k, i) / a(m, i)  
            a(k, :) = a(k, :) - f * a(m, :) 
            b(k)    = b(k)    - f * b(m)
          end do
        end do
        do i = size(b), 1, -1
          k = ipiv(i)  
          x(i) = ( b(k) - sum(a(k, i+1:) * x(i+1:)) ) / a(k, i)
        end do  
      end subroutine gauss0
    end module m_gauss

    program test1
      use m_timing 
      use m_gauss
      implicit none
      real, allocatable :: a(:, :), b(:), x(:)
      real, allocatable :: w(:, :), u(:)
      integer :: i, j, n = 2000
      call random_seed()      
      allocate( a(n, n), b(n), x(n) )
      call random_number(a)
      call random_number(b)
      w = a
      u = b
      call timimg('before gauss')
      call gauss(a, b, x)
      call timimg('after  gauss')
      print *, 'RMS of diff. ', sqrt(sum(abs(u - matmul(w, x))**2) / size(u))
      a = w
      b = u
      call timimg('before gauss')
      call gauss0(a, b, x)
      call timimg('after  gauss')
      print *, 'RMS of diff. ', sqrt(sum(abs(u - matmul(w, x))**2) / size(u))
    end program test1
```

##複素行列編
複素数 c = cr + i ci と、２x２実行列の同型性を利用して、実数用のルーチンで複素数計算をするテクニックがあります。複素変数での線型方程式 $Ax=b$ は、各変数の実部と虚部を r, i で表わすことにして、以下のような２倍サイズの実変数の問題に置き換えられます。
 
```math
\begin{pmatrix}
ar & -ai \\
ai & ar \\
\end{pmatrix}
\begin{pmatrix}
xr \\
xi \\
\end{pmatrix}
=
\begin{pmatrix}
br \\
bi \\
\end{pmatrix}
```

しかし、一般に線形方程式の解法は行列サイズ n の三乗に比例する計算量が必要なので、素朴にはこの場合の計算量は、元のサイズの実行列の場合にくらべて 2^3=8 倍の計算量（計算時間）になると予想されます。

一方で、複素数の積は、実部と虚部の積和計算になって素朴には、実数の積にくらべて
６倍程度の計算量になると予想されます。

これを実際の計算で確かめてみます。

###実行結果
2000 x 2000 の乱数で生成した複素行列で線型方程式を解いてみました。実際の計算は 4000 x 4000 サイズの実行列と、2000 x 2000 の複素行列で行っています。

計算結果を見ると、複素行列の場合、同サイズの実行列の時と比べて計算時間が３倍程度になっています。これは予想の６倍より少なくなっています。一方、２倍サイズの実行列の方は、半分サイズの実行列の場合に比べて２０倍以上の計算時間となっています。これは予想の８倍をはるかに超えています。

これらの原因は、まだ追及していないのでよく分かりませんが、２倍サイズの方は行列サイズが大きくなりすぎてキャッシュから溢れたとか、そういう事をとりあえず予想しています。複素行列が速い方は、Fortran では複素数は組み込み型なので、最適化が良くかかるのではないかと推測しています。

```
time =    0.00    0.00:before gauss
time =  277.09  277.09:after  gauss
 RMS of diff.   1.8752912E-04
time =    0.03  277.12:before gauss
time =   32.46  309.58:after  cgauss
 RMS of diff.   2.8650925E-04
```

###プログラム

```fortran
    module m_timing
      implicit none
    contains
      subroutine timimg(text)
        character(*), intent(in) :: text
        logical, save :: first = .true.
        integer, save :: ic0, ic1, ic2, icr
        if (first) then
          call system_clock(ic0, icr)
          ic1 = ic0
          first = .false.
        end if
        call system_clock(ic2, icr)
        print '(a, 2f8.2, 2a)', 'time =', real(ic2 - ic1) / icr, real(ic2 - ic0) / icr, ':', trim(text)
        ic1 = ic2
      end subroutine timimg
    end module m_timing
    
    module m_gauss
      implicit none
    contains
      subroutine gauss(a, b, x)
        real, intent(in out) :: a(:, :), b(:)
        real, intent(   out) :: x(:)
        integer :: i, j, k, ipiv(size(b))
        forall(i = 1:size(b)) ipiv(i) = i
        do i = 1, size(a, 2)
          k = maxloc(abs(a(ipiv(i:), i)), 1) + i - 1
          ipiv([i, k]) = ipiv([k, i]) ! swap pivot
          do j = i + 1, size(a, 1)     
            b(ipiv(j))    = b(ipiv(j))    - a(ipiv(j), i) / a(ipiv(i), i) * b(ipiv(i))
            a(ipiv(j), :) = a(ipiv(j), :) - a(ipiv(j), i) / a(ipiv(i), i) * a(ipiv(i), :) 
          end do
        end do
        do i = size(b), 1, -1
          x(i) = ( b(ipiv(i)) - sum(a(ipiv(i), i+1:) * x(i+1:)) ) / a(ipiv(i), i)
        end do  
      end subroutine gauss

      subroutine cgauss(a, b, x)
        complex, intent(in out) :: a(:, :), b(:)
        complex, intent(   out) :: x(:)
        integer :: i, j, k, ipiv(size(b))
        forall(i = 1:size(b)) ipiv(i) = i
        do i = 1, size(a, 2)
          k = maxloc(abs(a(ipiv(i:), i)), 1) + i - 1
          ipiv([i, k]) = ipiv([k, i]) ! swap pivot
          do j = i + 1, size(a, 1)     
            b(ipiv(j))    = b(ipiv(j))    - a(ipiv(j), i) / a(ipiv(i), i) * b(ipiv(i))
            a(ipiv(j), :) = a(ipiv(j), :) - a(ipiv(j), i) / a(ipiv(i), i) * a(ipiv(i), :) 
          end do
        end do
        do i = size(b), 1, -1 
          x(i) = ( b(ipiv(i)) - sum(a(ipiv(i), i+1:) * x(i+1:)) ) / a(ipiv(i), i)
        end do  
      end subroutine cgauss      
    
    end module m_gauss

    program test2
      use m_timing
      use m_gauss
      implicit none
      real, allocatable :: a(:, :), b(:), x(:), cr(:, :), ci(:, :), dr(:), di(:)
      real, allocatable :: w(:, :), u(:)
      
      complex, allocatable :: c(:, :), d(:), y(:)
      complex, allocatable :: z(:, :), v(:)
      
      integer :: i, j, n = 2000
      call random_seed()      
      allocate( a(2*n, 2*n), b(2*n), x(2*n), cr(n, n), ci(n, n), dr(n), di(n) )
      call random_number(cr)
      call random_number(ci)
      call random_number(dr)
      call random_number(di)
      a(   :n,    :n) =  cr 
      a(n+1: , n+1: ) =  cr
      a(   :n, n+1: ) = -ci
      a(n+1:2*n,  :n) =  ci
      b(   :n) = dr
      b(n+1: ) = di
      w = a
      u = b
      call timimg('before gauss')
      call gauss(a, b, x)
      call timimg('after  gauss')
      print *, 'RMS of diff. ', sqrt(sum(abs(u - matmul(w, x))**2) / size(u))
      
      allocate( c(n, n), d(n), y(n) )
      forall(i=1:n, j=1:n) c(i, j) = cmplx(cr(i, j), ci(i, j))
      forall(i=1:n       ) d(i)    = cmplx(dr(i)   , di(i)   )
      z = c
      v = d
      call timimg('before gauss')
      call cgauss(c, d, y)
      call timimg('after  cgauss')
      print *, 'RMS of diff. ', sqrt(sum(abs(v - matmul(z, y))**2) / size(v))
    end program test2
```