mvaxpy.f

来自「网络带宽测试工具」· F 代码 · 共 50 行

      Subroutine mvaxpy( a, b, c, lda, m, n )
! ----------------------------------------------------------------------
! --- Routine 'mvaxpy' does a matrix-vector multiplication 'Ab = c'
!     using a vector update implementation.
! ----------------------------------------------------------------------
      Use         numerics
      Implicit    None
      Integer  :: lda, m, n
      Real(l_) :: a(lda,n), b(n), c(m)

      Integer  :: i, chunk, iend, istart, j, l, me, nf, nthreads
      Integer  :: omp_get_num_threads, omp_get_thread_num
      Real(l_) :: lc(m)
! ----------------------------------------------------------------------
!$omp parallel do
      Do j = 1, m
         c(j)       = 0.0_l_
      End Do
      nf = n - Mod( n, 4 )
!$omp parallel private( chunk, me, nthreads, istart, iend, lc )
      nthreads = omp_get_num_threads()
      chunk    = m/nthreads
      me       = omp_get_thread_num()     
      istart   = me*chunk + 1
      If ( me < nthreads - 1 ) Then
         iend = ( me + 1 )*chunk
      Else
         iend = m
      End If
      Do i = 1, nf, 4
         Do j = istart, iend
            lc(j) =    a(j,i)  *b(i)   + a(j,i+1)*b(i+1)
     &               + a(j,i+2)*b(i+2) + a(j,i+3)*b(i+3)
         End Do
         Do j = istart, iend
            c(j) = c(j) + lc(j)
         End Do
      End Do
      Do i = nf + 1, n
         Do j = istart, iend
            lc(j) = a(j,i)*b(i)
         End Do
         Do j = istart, iend
            c(j) = c(j) + lc(j)
         End Do
      End Do
!$omp end parallel
! ----------------------------------------------------------------------
      End Subroutine mvaxpy