PRACE NPT

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

     &                   beta, y, incy, wrk )

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

! --- MVCPLN  performs one of the the following matrix-vector operations

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

!     where alpha and beta are scalars, x and y are vectors and A is an

!     m by n matrix.

!     The implementation is column oriented.

!     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  - The increment for the elements of X.  In this

!                     version, INCX must be 1 and is only here for

!                     compatibility with Level 2 BLAS routine DGEMV

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

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

!  INCY    - Input  - The increment for the elements of Y.  In this

!                     version, INCY must be 1 and is only here for

!                     compatibility with Level 2 BLAS routine DGEMV

!  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, 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( 'mvcpln', 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 j = 1, n

            temp = alpha * x(j)

            y(1:m) = y(1:m) + a(1:m,j)*temp

         End Do

      Else

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

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

        mybase = offset(me)

        Do j = 1, n

           temp = zero

           wrk(j) = zero

           Do i = 1, m

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

           End Do

           wrk(j) = wrk(j) + alpha*temp

        End Do

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

! ---   Combine partial results.

        type = MPI_Real8

        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 mvcpln