使用MKL中的有关函数计算共轭梯度法

program ConjugateGradientMethods
        use blas95
        implicit none
        integer :: i
        integer, parameter :: n = 5
        real(kind=8),allocatable :: A(:,:)
        real(kind=8),allocatable :: b(:)
        real(kind=8) :: alpha = 0.d0, beta = 0.d0
        real(kind=8) :: x(n), d(n), r(n), tmp
  
        !// 默认初始值x0 = 0; 则 d = b - A * x0 = b
        x = 0.d0; d = 0.d0; r = 0.d0
        allocate( A(n,n), b(n) )
        
        open( 101, file = 'A.txt' )
        read( 101,* ) A
        read( 101,* ) b
        close( 101 )
        
        !$-----------------------------------------------$
        !|                   A x = b                     |
        !$-----------------------------------------------$
        !| 0.7719   -0.0133    0.0043   -0.1540   -0.0008|
        !|-0.0133    0.3392    0.1065   -0.0090   -0.0210|
        !| 0.0043    0.1065    0.2976    0.0029   -0.0568|
        !|-0.1540   -0.0090    0.0029    0.8960   -0.0006|
        !|-0.0008   -0.0210   -0.0568   -0.0006    0.3103|
        !|-----------------------------------------------$
        !| 0.1380    0.8440    0.8375    3.4180    1.3357|
        !$-----------------------------------------------$
        
        d = b
        r = d

        i = 1
        do 
                if ( minval(abs(r)) < 1d-10 .or. i > 100 ) exit

                tmp = dot(r, r)

                call gemv( A, d, b )
                alpha = tmp / dot( d, b )
                call axpy( d, x, alpha )  !// x = x + alpha * d  <update x>
                call axpy( b, r, -alpha ) !// r = r - alpha * b  <update r>

                beta = dot(r, r) / tmp
                call scal( d, beta )
                call axpy( r, d, 1.d0 )   !// d = beta * d + r   <update d>
    
                i = i + 1
        end do
        deallocate( A, b )
  
        do i = 1, n
                write ( *,* ) x(i)
        end do
end program ConjugateGradientMethods