PRACE NPT

      Program mod2a

! **********************************************************************

! *** This program is a distributed-memory version of                ***

! *** EuroBen Benchmark program MOD2A (generalized matrix-vector     ***

! *** multiplication).                                               ***

! ***                                                                ***

! *** Copyright: European Benchmark Group p/o                        ***

! ***            Computational Physics Dept. Utrecht                 ***

! ***            P.O. Box 80125                                      ***

! ***            3508 TD Utrecht                                     ***

! ***            The Netherlands                                     ***

! ***                                                                ***

! *** Author of this program: Aad J. van der Steen                   ***

! *** Date                    04/29/1998                             ***

! *** Based on the shared-memory version of Ruud van der Pas,        ***

! ***                                           11/24/1988           ***

! *** and Peter de Rijk                       , 02/01/1993           ***

! *** Improved communication                  , 05/14/1996           ***

! *** Total rework                            , 04/29/1998           ***

! **********************************************************************

! --- Version 3.3

!

! --- Purpose of program Mod2a

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

!     This program measures the performance of the matrix-vector product

!     y = A x  and  y = A' x, for a full  n x n  matrix A  with a vector

!     of length  n.

!     Several variants are examined.

!

! --- Remarks

!     -------

!     It is possible to measure the performance of a special library

!     routine for matrix vector multiplication.

!     If such a routine is available, the recipe is:

!       - Replace the string DGEMV in character string LNAME by the

!         appropriate name.

!       - Set parameter ILIB to a non-zero integer value.

!       - Activate in subroutine LIBVER the call to DGEMV or replace

!         this statement with the appropriate call to the library

!         routine.

!     Warning: To prevent making mistakes, there is NO source version of

!              DGEMV included in this program.

!              The default situation is ILIB=0 i.e. no library version

!              present!

!

! --- A EuroBen encapsulation routine is used for setting the network 

!     (csetup). In addition, the EuroBen status routine 'state', and

!     the routine for accessing. Routine 'state' should be customized

!     for the machine at hand.

!

! --- Presently, the maximum number of nodes allowed is 2048. To

!     increase this bound, increase the value of parameter 'maxnod' in 

!     in all relevant routines.

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

      Use                   dist_module

      Use                   numerics

      Integer, Parameter :: nmax = 2500 

      Integer, Parameter :: ncases = 4, lda = nmax, mxcase = 50,

     &                      ilib = 0

      Character          :: modnam*8, text(ncases+1)*20, lname*43

      Integer            :: icase, info, ndim(mxcase)

      Integer            :: i, j, k, n, m, mbase, nrep, ipoint,

     &                      lchk(nmax)

      Real(l_)           :: tiperf((ncases+1)*mxcase,4), extime, perfor

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

     &                      ychk(nmax)

      Real(l_)     :: start_time, end_time

      External              mvrpln, mvrur4, mvcpln, mvcur4

      Data modnam /        'mod2a   ' /

!                           123456789 123456789 1

      Data text   /        'Rows,    plainly    ' ,

     &                     'Rows,    unrolled 4 ' ,

     &                     'Columns, plainly    ' ,

     &                     'Columns, unrolled 4 ' ,

     &                     'Library routine(s)  ' /

!

!                           123456789 123456789 123456789 123456789 1234

      Data lname  /        'Library routine(s) used: DGEMV           ' /

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

! --- Set up communication network and print test status.

      Call csetup

      If ( me == 0 ) Call state (modnam)

      start_time = MPI_Wtime()

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

! --- Open file containing the repetition factor and the problem size

!     for each problem. Read and compute until EOF (maximum no. of

!     cases allowed is mxcase = 50). The maximum problem size

!     nmax = 2500.

      icase = 0

      Open( 1, File = 'mod2a.in' )

   10 Read( 1, *, End = 20 ) nrep, n

      icase = icase + 1

      If ( icase > mxcase ) Then

         If ( me == 0 ) Print *, 'More than mxcase = ', mxcase, 

     &                    ': increase mxcase'

         Go To 20

      End If

      ndim(icase) = n

      If ( n > nmax ) Then

         If ( me == 0 ) Print *, 'n > nmax = ', nmax,

     &                           ': increase nmax'

         Go To 20

      End If

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

! --- Distribute a as evenly as possible over the available processors.

      Call evdist( n )

      Call bsaddr

      m = sizes(me)

! --- Generate the matrix A and vector x.

      Call matgen( a, x, lda, n )

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

! --- For checking the correctness we generate the solution vector by

!     independent means:

      Call mkbnds( a, lda, ychk )

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

!                                                ** Row oriented variant

      ipoint = icase

      Call timing( 'N', m, n, a, lda, x, y, mvrpln,

     &             nrep, extime, perfor, wrk )

! --- Check correctness:

       Call check( 'mvrpln', 'N', y, ychk, lchk )

       tiperf(ipoint,1) = extime

       tiperf(ipoint,2) = perfor

       Call timing( 'T', m, n, a, lda, x, y, mvrpln,

     &              nrep, extime, perfor, wrk )

! --- Check correctness:

      Call check( 'mvrpln', 'T', y, ychk, lchk )

!

      tiperf(ipoint,3) = extime

      tiperf(ipoint,4) = perfor

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

!                            ** Row oriented variant with loop unrolling 

      ipoint = icase + mxcase

      Call timing( 'N', m, n, a, lda, x, y, mvrur4,

     &             nrep, extime, perfor, wrk )

! --- Check correctness:

      Call check( 'mvrur4', 'N', y, ychk, lchk )

      tiperf(ipoint,1) = extime

      tiperf(ipoint,2) = perfor

      Call timing( 'T', m, n, a, lda, x, y, mvrur4,

     &             nrep, extime, perfor, wrk )

! --- Check correctness:

      Call check( 'mvrur4', 'T', y, ychk, lchk )

      tiperf(ipoint,3) = extime

      tiperf(ipoint,4) = perfor

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

!                                             ** Column oriented variant

      ipoint = icase + 2*mxcase

      Call timing( 'N', m, n, a, lda, x, y, mvcpln,

     &             nrep, extime, perfor, wrk )

! --- Check correctness:

      Call check( 'mvcpln', 'N', y, ychk, lchk )

      tiperf(ipoint,1) = extime

      tiperf(ipoint,2) = perfor

      Call timing( 'T', m, n, a, lda, x, y, mvcpln,

     &             nrep, extime, perfor, wrk )

! --- Check correctness:

      Call check ( 'mvcpln', 'T', y, ychk, lchk )

      tiperf(ipoint,3) = extime

      tiperf(ipoint,4) = perfor

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

!                         ** Column oriented variant with loop unrolling

      ipoint = icase + 3*mxcase

      Call timing( 'N', m, n, a, lda, x, y, mvcur4,

     &             nrep, extime, perfor, wrk )

! --- Check correctness:

      Call check( 'mvcur4', 'N', y, ychk, lchk )

      tiperf(ipoint,1) = extime

      tiperf(ipoint,2) = perfor

      Call timing( 'T', m, n, a, lda, x, y, mvcur4,

     &             nrep, extime, perfor, wrk )

! --- Check correctness:

      Call check( 'mvcur4', 'T', y, ychk, lchk )

      tiperf(ipoint,3) = extime

      tiperf(ipoint,4) = perfor

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

!                                             ** Special library variant

      If ( ilib /= 0 ) Then

         ipoint = icase + 4*mxcase

! --- Special library version of matrix vector multiplication

         Call libver( 'N', m, n, a, lda, x, y, 

     &                nrep, extime, perfor, wrk )

! --- Check correctness:

         Call check( 'libver', 'N', y, yck, lchk )

         tiperf(ipoint,1) = extime

         tiperf(ipoint,2) = perfor

         Call libver( 'T', m, n, a, lda, x, y, 

     &                nrep, extime, perfor, wrk )

! --- Check correctness:

         Call check( 'libver', 'T', y, yck, lchk )

         tiperf(ipoint,3) = extime

         tiperf(ipoint,4) = perfor

      EndIf

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

! >>> Get new case at label 10 <<<

      Go To 10

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

!                                                  ** Print the results.

   20 If ( me == 0 ) Then

         Print 9010, nodes, nrep

         Do k = 1,ncases

            ipoint = (k-1)*mxcase + 1

            Print 9020, text(k), ndim(1), ndim(1), 

     &                  (tiperf(ipoint,j), j=1,4)

            Do i = 2, icase

               ipoint = ipoint + 1

               Print 9030, ' ',ndim(i), ndim(i), (tiperf(ipoint,j),

     &                     j=1,4)

            End Do

         End Do

!

         If ( ilib /= 0 ) Then

            ipoint = ncases*mxcase + 1

            Print 9020, text(ncases+1), ndim(1), ndim(1), 

     &                 (tiperf(ipoint,j), j=1,4)

            Do i = 2,nt

               ipoint = ipoint + 1

               Print 9030, ' ', ndim(i), ndim(i), (tiperf(ipoint,j),

     &                     j=1,4)

            End Do

            Print 9040, lname

         End If

      End If

      end_time =  MPI_Wtime() - start_time

      If (me == 0) Then

         Write(6,22) 'Walltime: ', end_time, " s"

 22      Format(A,F9.3,A)

      End If

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

!                                               ** Exit network orderly.

      Call MPI_Finalize( info )

      If ( me == 0 ) Then

      Write(6,*) 'Program terminated normally'

      End If

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

!                                                            ** Formats.

 9010 Format( //79('-')/

     &'Generalized matrix-vector multiplication'/

     &4X,'y := alpha*A*x + beta*y   or',

     &3X,'y := alpha*A''*x + beta*y',/

     &'Several variants are measured, No. of procs. = ', i3 /

     &79('-')//

     &'Initial repetition factor used in timings: ',I3//

     &'Variant, i.e. array',18X,'No Transpose',12X,'Transpose'/

     &'A is accessed by',4X,4X,'M',4X,'N',

     &    2(4X,'seconds',4X,'Mflop/s ')/

     &79('-') )

 9020 Format(/A20,2(1X,I4),1P,2(3X,E9.3,2X,E9.3) )

 9030 Format( A20,2(1X,I4),1P,2(3X,E9.3,2X,E9.3) )

 9040 Format( A43 )

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

      End Program mod2a