*
* ACML version 3.5 Copyright AMD,NAG 2006
*
      PROGRAM MAIN
      IMPLICIT         NONE
*     .. Parameters ..
      INTEGER          MAXN
      PARAMETER        (MAXN=5000)
      INTEGER          LDA, NN
      PARAMETER        (LDA=MAXN,NN=11)
      DOUBLE PRECISION ALIM, BLIM
      PARAMETER        (ALIM=1.0D0,BLIM=2.0D0)
*     .. Local Scalars ..
      INTEGER          J, M, N
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,MAXN)
      INTEGER          ENS(NN), IPIVA(LDA)
*     .. External Subroutines ..
      EXTERNAL         DOTIME
*     .. Data statements ..
      DATA             ENS/100, 400, 800, 1200, 1600, 2000, 2400, 2800,
     *                 3200, 3600, 4000/
*     .. Executable Statements ..
      CONTINUE
*
      DO 20 J = 1, NN
*
         M = ENS(J)
*
*        To time different data, change ENS, NN, and MAXN, above,
*        accordingly.
*
         N = M
*
         IF (M.GT.MAXN) THEN
            WRITE (*,FMT=*) 'M too large; increase MAXN.'
         ELSE
            CALL DOTIME(M,N,A,LDA,IPIVA,ALIM,BLIM)
         END IF
*
   20 CONTINUE
*
      END
*
*-----------------------------------------------------------------------
*-----------------------------------------------------------------------
*
      SUBROUTINE DOTIME(M,N,A,LDA,IPIVA,ALIM,BLIM)
      IMPLICIT         NONE
*     .. Scalar Arguments ..
      DOUBLE PRECISION ALIM, BLIM
      INTEGER          LDA, M, N
*     .. Array Arguments ..
      DOUBLE PRECISION A(LDA,*)
      INTEGER          IPIVA(*)
*     .. Local Scalars ..
      DOUBLE PRECISION DNFLOP, MFX, T0X, T1X, TX
      INTEGER          INFO
*     .. External Functions ..
      DOUBLE PRECISION DSECND
      EXTERNAL         DSECND
*     .. External Subroutines ..
      EXTERNAL         DGETRF, RANDMAT
*     .. Intrinsic Functions ..
      INTRINSIC        DBLE
*     .. Executable Statements ..
      CONTINUE
*
*     Number of MFLOPs in a DGETRF of size M-by-N:
*
      IF (N.LT.M) THEN
         DNFLOP = (DBLE(M)*(DBLE(N)**2)-(DBLE(N)**3)/3.0D0)/1.0D6
      ELSE IF (N.EQ.M) THEN
         DNFLOP = ((2.0D0/3.0D0)*DBLE(N)**3)/1.0D6
      ELSE
         DNFLOP = (DBLE(N)*(DBLE(M)**2)-(DBLE(M)**3)/3.0D0)/1.0D6
      END IF
*
      CALL RANDMAT(ALIM,BLIM,M,N,A,LDA)
*
*     Using DGETRF
*
      T0X = DSECND()
      CALL DGETRF(M,N,A,LDA,IPIVA,INFO)
      T1X = DSECND()
      TX = T1X - T0X
*
*     Print the results
*
      IF (INFO.EQ.0) THEN
         IF (TX.LE.0.0D0) THEN
            MFX = 0.0D0
         ELSE
            MFX = DNFLOP/TX
         END IF
         WRITE (*,FMT=*) N, MFX
      ELSE
         WRITE (*,FMT=*) 'DGETRF failed with INFO = ', INFO
      END IF
*
      RETURN
      END
*
      SUBROUTINE RANDMAT(ALIM,BLIM,M,N,A,LDA)
*     .. Parameters ..
      INTEGER          MSEED, MSTATE
      PARAMETER        (MSEED=624,MSTATE=633)
*     .. Scalar Arguments ..
      DOUBLE PRECISION ALIM, BLIM
      INTEGER          LDA, M, N
*     .. Array Arguments ..
      DOUBLE PRECISION A(LDA,N)
*     .. Local Scalars ..
      INTEGER          GENID, I, INFO, LSEED, LSTATE, SUBID
*     .. Local Arrays ..
      INTEGER          SEED(MSEED), STATE(MSTATE)
*     .. External Subroutines ..
      EXTERNAL         DRANDINITIALIZE, DRANDUNIFORM
*     .. Executable Statements ..
      CONTINUE
*
*     Use the Mersenne twister generator as the base generator
*
      GENID = 3
      SUBID = 1
*
      LSTATE = MSTATE
      LSEED = 1
      SEED(1) = 1071958
*
*     Initialize the base generator
*
      CALL DRANDINITIALIZE(GENID,SUBID,SEED,LSEED,STATE,LSTATE,INFO)
*
      DO 20 I = 1, N
*
*        Generate a matrix from a uniform U(ALIM,BLIM) distribution
*
         CALL DRANDUNIFORM(M,ALIM,BLIM,STATE,A(1,I),INFO)
*
   20 CONTINUE
*
      RETURN
      END
