Householder法による三重対角化

昔Householder法とQR法のサブルーチンを外人からもらいましたが、EISPACKの tred と tql で F66 で書かれていました。のちに Numerical Recipes の F77 版に置き換えたけれども、基本的に LAPACK のものだと思います。微妙に変えてあって、tql は旧版で収束しない時でも計算できるようになったけれど収束判定が厳しくなりすぎて遅くなった記憶が。

tred は昔のプログラムなので、記憶領域を使わないように色々芸当を使っていますが、策士策に溺れるような感じがしなくもないような・・ワーク配列を取らないで、元の行列の空きスペースを利用するようにしています。７７時代はワーク配列を呼び出し側で用意しなければならないので、便利ではあるんですが、おかげで解読しにくくなっています。ここでは、ワーク配列をベクトル２本分用意して、できるだけ数式通りになる様にしました。配列はメモリーが連続アクセスされるようにしています。また、変換行列を求める部分を別ルーチンにしました。

実行結果

Householder法で対称行列を三重対角化して、対角要素と非対角要素を出力します。また元の行列と三重対角化された行列を結ぶ変換行列を求め、実際に変換してみて、三重対角行列が得られることを確かめます。さらに対角要素の総和を元の行列、三重対角行列で求めて比較しています。

ソース・プログラム

module m_eigen
    implicit none
    integer, parameter :: kd = kind(0.0d0)
 contains
 
    subroutine hh(a, e, d)
      real(kd), intent(in out):: a(:, :)
      real(kd), intent(out) :: e(:), d(:)
      real(kd) :: p(size(a, 1)), q(size(a, 1))
      real(kd) :: s, t, ak, alpha, scale  
      integer :: i, k, k1, n
      n = size(a, 2)
      do k = 1, n - 2
        k1 = k + 1
        scale = maxval(abs(a(k1:, k)))
        a(k1:, k) = a(k1:, k) / scale
        t = dot_product(a(k1:, k), a(k1:, k))
        ak = a(k1, k)
        s = sign(sqrt(t), ak)
        alpha = 1.0_kd / (t + ak * s) 
        e(k) = a(k, k)  
        d(k) = -s * scale 
        a(k, k)  = alpha   ! data for Eigen vectors
        a(k1, k) = ak + s  ! .. w = a - b,  alpha = 1 / |w|
        forall(i = k1:n) p(i) = alpha * dot_product( a(i, k1:), a(k1:, k) ) ! p = alpha*A*w
        t = 0.5_kd * alpha * dot_product( p(k1:), a(k1:, k) )
        q(k1:) = p(k1:) - t * a(k1:, k)                         ! q = p - 1/2alpha*w*p^t*w
        forall(i = k1:n) a(k1:, i) = a(k1:, i) - a(k1:, k) * q(i) - q(k1:) * a(i, k) 
      end do
      e(n - 1) = a(n - 1, n - 1)
      e(n)     = a(n, n)
      d(n - 1) = a(n, n - 1) 
      a(n - 1, n - 1) = 1.0_kd ! for Eigen vectors
      a(n, n) = 1.0_kd         !..
      a(n - 1, n) = 0.0_kd     !..
      a(n, n - 1) = 0.0_kd     !..
      return    
    end subroutine hh
    
    subroutine qp(a)
      real(kd), intent(inout) :: a(:, :)
      real(kd) :: g
      integer :: i, j, n
      n = size(a, 1)
      do i = n - 2, 1, -1
        do j = i + 1, n                                       ! A = (I - alpha*u*u^t) * A
          g = a(i, i) * dot_product( a(i + 1:, i), a(i + 1:, j) ) ! alpha * u^t * A
          a(i + 1:, j) = a(i + 1:, j) - a(i + 1:, i) * g         ! A = A - u * g  
        end do
        a(i, i) = 1.0_kd
        a(i, i + 1:) = 0.0_kd
        a(i + 1:, i) = 0.0_kd
      end do
      return
    end subroutine qp
  end module m_eigen

  program Householder
    use m_eigen 
    implicit none
    integer, parameter :: n = 15
    real(kd) :: a(n, n), c(n, n), cc(n, n), e(n), d(n)
    integer :: i
    call random_seed()
    call random_number(a)

    a = - ( a + transpose(a) ) ! symmetrize
    c = a
    call hh(a, e, d)
    call qp(a)
    cc = matmul(transpose(a), matmul(c, a)) 

    do i = 1, n
      print '(2f6.2)', e(i), d(i)
    end do
    print *
    do i = 1, n
      print '(100f6.2)', cc(i, :)
    end do

    print *, ' trace check ' 
    print *, sum( [(c(i, i), i = 1, n)] )
    print *, sum( [(cc(i, i), i = 1, n)] )
    print *, sum( e )
    stop
  end program Householder