STAMP CPUTIME

[Fortran95]メモ帳

 module m_stamp
    implicit none
  contains
    subroutine stamp(text)
      character (len = *), intent(in), optional :: text
      character (len = 40) :: buff
      character (len = 10) :: date, time
      real, save :: time0, time1
      logical, save :: qfirst = .true.
      integer, save :: icount0, icount1, irate
      if (qfirst) then
        call cpu_time(time0)
        call system_clock(icount0, irate)
        qfirst = .false.
      end if
      if (present(text)) then
        buff = ': ' // trim(text)
      else
        buff = ''
      end if
      call date_and_time( date, time )
      call cpu_time(time1)
      call system_clock(icount1)
      print '(6a, 1x , 6a, f8.1, 2a, f8.1, 2a)', '#', &
            date(1:4), '-', date(5:6), '-', date(7:8), time(1:2), ':', time(3:4), ':', time(5:6), &
            ' cpu time', time1 - time0, 'sec ', 'elapsed time ', real(icount1 - icount0) / irate, 'sec ', buff 
      return
    end subroutine stamp
  end module m_stamp

  module m_cg
    implicit none
    integer, parameter :: kd = selected_real_kind(12) !double precision
  contains
    subroutine cg(a, b, x) ! Conjugate Gradient method
      real(kd), intent(in) :: a(:, :), 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 - matmul(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 = matmul(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_cg
  
  program stamps
    use m_stamp
    use m_cg
    implicit none
    integer, parameter :: ns = 2000 
    real(kd) :: a(ns, ns), b(ns), x(ns)
    integer :: i, j
  ! Make matrix 
    call stamp()
    do i = 1, ns
      b(i) = 1.0_kd
      do j = i, ns ! Givens Matrix,  x = (1, 0, 0, ... , 0)t
        a(i, j) = real( 2 * MIN(i, j) - 1, kd )  
        a(j, i) = a(i, j)
      end do
    end do
    call stamp('before cg')

    ! SOLVE Ax = b
    call random_seed() 
    call random_number(x) ! starting vector
    call cg(a, b, x)     
    call stamp('after  cg')

    x = matmul(a, x) - b
    print *, 'difference (RMS) =', sqrt( dot_product(x, x) )
    call stamp()
    stop
  end program stamps