【メモ帳】abstract type の中に concrete な routine をぶら下げる

記事を移動し独立させました

昔のプログラムが、今の ifort を通らないので修正。

元々 abstract type に concrete なルーチンをくっつけていて違和感があったので、書き直してみました。

どう直せばいいのか少し苦労しましたが、メンバー要素として type bound procedure pointer にすることで解決しました。

fortran66.hatenablog.com

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)

    module m_base
        implicit none
        private
        integer, parameter, public :: kd = kind(1.0d0)
!
        type, abstract, public :: t_matrix
            procedure(cg_sub), pointer :: cg => cg_sub
        contains
            procedure(prd_fun), deferred :: product_abst
            generic, public :: operator(*) => product_abst
        end type t_matrix
    
        abstract interface
            function prd_fun(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 prd_fun
        end interface     

    contains

        subroutine cg_sub(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
        end subroutine cg_sub
    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, public :: product_abst => 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)
        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 :: product_abst => 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 
        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)
    end program CG