密行列と粗行列のCG法を単一のルーチンで記述

密行列と粗行列のCG法（共役勾配法）による線型方程式の解法を、オブジェクト指向機能を用いて単一のルーチンで記述してみました。行列とベクトルの積をそれぞれの行列形式について定義し、それを「＊」演算子に総称化させることで可能になります。

abstract type は単なるスケルトン記述しか許されないのかと思っていましたが、具体的なルーチンも記述可能なようです。

実行結果

 dense  matrix |a * x - b| = 4.268173152696667E-007
 sparse matrix |q * y - b| = 4.268173152696667E-007
 x==y? T
続行するには何かキーを押してください . . .

ソース・プログラム

    module m_base
      implicit none
      private
      integer, parameter, public :: kd = kind(1.0d0)
!
      type, abstract, public :: t_matrix
      contains
        procedure, public :: cg
        procedure(op2), deferred, private :: prdct
        generic, public :: operator(*) => prdct
      end type t_matrix
      
      abstract interface
        function op2(mat, vec) result(res) 
          import :: t_matrix, kd 
          class(t_matrix), intent(in) :: mat
          real(kd), intent(in) :: vec(:) 
          real(kd) :: res(size(vec))
        end function op2
      end interface     

    contains

      subroutine cg(A, b, x) ! Conjugate Gradient method
        class(t_matrix), intent(in) :: A
        real(kd), intent(in) :: b(:)
        real(kd), intent(in out) :: x(:)
        real(kd) :: alpha, beta, r2n, r2d, eb2, p(size(b)), q(size(b)), r(size(b))
        integer :: i
        r = b - A * x
        p = r
        r2n = dot_product(r, r)
        eb2 = epsilon(b) * dot_product(b, b)
        do i = 1, size(b)
          if ( r2n < eb2 ) exit
          r2d = r2n
          q = A * p
          alpha = dot_product(p, r) / dot_product(p, q)
          x = x + alpha * p
          r = r - alpha * q
          r2n = dot_product(r, r)
          beta = r2n / r2d
          p = r + beta * p            
        end do
        return
      end subroutine cg
    end module m_base

    module m_dense
      use m_base
      implicit none
      private
      type, extends(t_matrix), public :: t_dense_matrix
        real(kd), allocatable :: x(:, :)
      contains
        procedure :: prdct => product_dense  
      end type t_dense_matrix

    contains
    
      function product_dense(mat, vec) result(res)
        class(t_dense_matrix), intent(in) :: mat
        real(kd), intent(in) :: vec(:)
        real(kd) :: res(size(vec))
        res = matmul(mat%x, vec)
        return
      end function product_dense
    end module m_dense
    
    module m_sparse
      use m_base
      implicit none
      private
      type :: t_sparse_vec
        integer , allocatable :: idx(:)
        real(kd), allocatable :: x(:)
      end type t_sparse_vec
      
      type, extends(t_matrix), public :: t_sparse_matrix
        type(t_sparse_vec), allocatable :: row(:) 
      contains
        procedure :: prdct => product_sparse  
      end type t_sparse_matrix

    contains
      
      function product_sparse(mat, vec) result(res)
        class(t_sparse_matrix), intent(in) :: mat
        real(kd), intent(in) :: vec(:)
        real(kd) :: res(size(vec))
        real(kd) :: tmp
        integer :: i, k
        do i = 1, size(mat%row)
          tmp = 0.0_kd
          do k = 1, size(mat%row(i)%idx)
            tmp = tmp + mat%row(i)%x(k) * vec(mat%row(i)%idx(k))
          end do    
          res(i) = tmp
        end do 
        return
      end function product_sparse
    end module m_sparse
    
    program CG
      use m_base
      use m_dense
      use m_sparse
      implicit none
      type(t_dense_matrix ) :: a
      type(t_sparse_matrix) :: q
      real(kd), allocatable :: x(:), y(:), b(:)
      integer :: i, j, n = 1000  
      !
      allocate( a%x(n, n), b(n), x(n) ) ! dense
      allocate( q%row(n)       , y(n) ) ! sparse
      !
      do i = 1, n
        b(i) = 1.0_kd
        do j = i, n ! x = (1, 0, 0, ... , 0)t
          a%x(i, j) = real( 2 * MIN(i, j) - 1, kd )  
          a%x(j, i) = a%x(i, j)
        end do
      end do      
      do i = 1, n                             ! copy dense-mat to sparse-mat 
        q%row(i)%idx = [(j, j = 1, n)]  
        q%row(i)%x   = a%x(i, :)
      end do      
      !
      call random_seed()
      call random_number(x) ! starting vector
      y = x
      call a%cg(b, x)
      call q%cg(b, y)
      print *, 'dense  matrix |a * x - b| =', norm2(a * x - b) ! norm2 : f2008 
      print *, 'sparse matrix |q * y - b| =', norm2(q * y - b)
      print *, 'x==y?', all(x==y)
      stop
    end program CG