PRACE NPT

      Subroutine mvrpln( trans, m, n, alpha, a, lda, x, incx,

     &                   beta, y, incy, wrk )

! ----------------------------------------------------------------------

! --- Row oriented implementation of matrix-vector multiplication 

!     y = alpha*A*x + beta*y   or   y = alpha*A'*x + beta*y

!     The calling sequence is based on the Level 2 BLAS routine DGEMV

!  Parameters

!  ----------

!  TRANS   - Input  - Specifies the operation to be performed:

!                     TRANS = 'N' implies  y = alpha*A*x + beta*y.

!                     TRANS = 'T' implies  y = alpha*A'*x + beta*y.

!  M, N    - Input  - Dimension of the problem.

!  ALPHA   - Input  - Specifies the scalar alpha.

!  A(M,N)  - Input  - Contains the matrix A.

!  LDA     - Input  - Leading Dimension of A.

!  X(*)    - Input  - Contains the vector x.

!  INCX    - Input  - Specifies the increment for the elements of X.

!                     In this implementation, only INCX = 0 is allowed.

!                     This parameter is here only for compatibility 

!                     with Level 2 BLAS routine DGEMV only.

!  BETA    - Input  - Scalar beta; BETA = 0, then Y need not be set.

!  Y(*)    - Output - Contains result vector y

!

!  INCY    - Input  - Specifies the increment for the elements of Y.

!                     In this implementation, only INCY = 0 is allowed.

!                     This parameter is here only for compatibility 

!                     with Level 2 BLAS routine DGEMV only.

!  WRK(*)  - Workspace. NOT compatible with DGEMV!

! ----------------------------------------------------------------------

      Use                    dist_module

      Use                    numerics

      Character*1         :: trans

      Integer             :: m, n, lda, incx, incy

      Real(l_)            :: alpha, beta

      Real(l_)            :: a(lda,*), x(*), y(*), wrk(*)

! ----------------------------------------------------------------------

! --- Communication types:

      Include                'mpif.h'

      Integer             :: com, type

! --- Local constants

      Real(l_), Parameter :: zero = 0.0_l_, one = 1.0_l_

! --- Local variables

      Integer             :: i, info, j, jt, k, kbase, mybase

      Real(l_)            :: temp

! ----------------------------------------------------------------------

! --- Test the input parameters

      info = 0

      If ( trans /= 'N' .AND. trans /= 'T' ) Then

         info = 1

      Else If ( m < 0 ) Then

         info = 2

      Else If ( n < 0 ) Then

         info = 3

      Else If( lda < MAX ( 1, m ) ) Then

         info = 6

      Else If ( incx /= 1 ) Then

         info = 8

      Else If ( incy /= 1 ) Then

         info = 11

      End If

      If ( info /= 0 ) Then

         Call eberr ( 'mvrpln', info )

         Return

      End If

! ----------------------------------------------------------------------

! --- Start the operations; First form y := beta*y

      If ( beta /= one ) Then

         If ( beta == zero ) Then

            y(1:m) = zero

         Else

            y(1:m) = beta*y(1:m)

         End If

      End If

      If ( alpha == zero ) Return

      If ( trans == 'N' ) Then

! ----------------------------------------------------------------------

! --- Form  y := y + alpha*A*x

         Do i = 1, m

            temp = zero

            Do j = 1, n

               temp = temp + a(i,j) * x(j)

            End Do

            y(i) = y(i) + alpha*temp

         End Do

      Else

! ----------------------------------------------------------------------

! --- Form  y := y + alpha*A'*x; generate partial results.

        mybase = offset(me)

        wrk(1:n) = zero

        Do i = 1, m

           temp = alpha*x(mybase+i)

           wrk(1:n) = wrk(1:n) + a(i,1:n)*temp

        End Do

! ----------------------------------------------------------------------

! --- Combine partial results.

        type = MPI_Double_Precision

        com  = MPI_Comm_World

        Do k = 0, nodes - 1

           kbase = offset(k) + 1

           Call MPI_Reduce( wrk(kbase), y, sizes(k), type, MPI_Sum, k, 

     &                      com, info ) 

        End Do

      Endif

! ----------------------------------------------------------------------

      End Subroutine mvrpln