# To unbundle, sh this file
echo check.f 1>&2
cat >check.f <<'End of check.f'
      subroutine check(ta,lda,ldb,ldc,n1,n2,n3,a,b,c,cc)
*
*     Jack Dongarra and Dan Sorensen, 1/26/87.
*     Argonne National Laboratory
*
c
      logical ta
      integer lda,ldb,ldc,n1,n2,n3
      double precision a(lda,*),b(ldb,*),c(ldc,*),cc(ldc,*)
      double precision cnorm,rnorm,maxval
c
      maxval = 1.0e-10
c
      if( ta )then
         call zero(ldc,n2,n3,cc)
         n = n2
         do 30 j = 1,n3
            do 20 i = 1,n2
               do 10 k = 1,n1
                  cc(i,j) = cc(i,j) + a(k,i)*b(k,j)
   10          continue
   20       continue
   30    continue
      else
         call zero(ldc,n1,n3,cc)
         n = n1
         do 60 j = 1,n3
            do 50 k = 1,n2
               do 40 i = 1,n1
                  cc(i,j) = cc(i,j) + a(i,k)*b(k,j)
   40          continue
   50       continue
   60    continue
      end if
      cnorm = 0.0
      rnorm = 0.0
      do 80 j = 1,n3
         do 70 i = 1,n
            rnorm = max(rnorm,abs(c(i,j) - cc(i,j)))
            cnorm = max(cnorm,abs(cc(i,j)))
   70    continue
   80 continue
      if( rnorm.gt.maxval*cnorm ) then
         write(6,*)' error in matrix multiply'
         write(6,*)' norm of the difference is = ',rnorm
      endif
      return
      end
End of check.f
echo dgemm.f 1>&2
cat >dgemm.f <<'End of dgemm.f'
      subroutine dgemm ( transa, transb, m, n, p, alpha, a, lda, b, ldb,
     $                   beta, c, ldc )
*
*     Jack Dongarra and Dan Sorensen, 1/26/87.
*     Argonne National Laboratory
*
      character*1        transa, transb
      integer            m, n, p, lda, ldb, ldc
      double precision   alpha, beta
      double precision   a( lda, * ), b( ldb, * ), c( ldc, * )
*
*     dgemm   performs one of the matrix-matrix operations given by
*
*        C := alpha*A*B + beta*C,     transa = 'n', transb = 'n',
*
*           where  A is m by n,  B is n by p  and  C is m by p,
*
*        C := alpha*A'*B + beta*C,    transa = 't', transb = 'n',
*
*           where  A is m by n,  B is m by p  and  C is n by p,
*
*        C := alpha*A*B' + beta*C,    transa = 'n', transb = 't',
*
*           where  A is m by n,  B is p by n  and  C is m by p,
*
*        C := alpha*A'*B' + beta*C,   transa = 't', transb = 't'.
*
*           where  A is m by n,  B is p by m  and  C is n by p.
*
*     For both  transa  and  transb  'c'  has the same meaning as  't'.
*
*
*
      intrinsic          matmul, transpose
      logical            ta    , tb
*
*     if lda ldb and ldc different a problem
      ta = ( transa.eq.'t' ).or.( transa.eq.'c' )
      tb = ( transb.eq.'t' ).or.( transb.eq.'c' )
      if( transb.eq.'n' )then
         if( ( transa.eq.'n' ).and.( alpha.eq.1.0d+0 ).and.
     $       ( beta.eq.1.0d+0 )                             )then
            lvect = 64
            icache = 10000
            lblock = icache/lvect
            do 50 ii = 1,m,lvect
               iiend = min(ii+lvect-1,m)
               do 40 kk = 1,n,lblock
                  kkend = min(kk+lblock-1,n)
cvd$ cncall
                  do 30 j = 1,p,4
                     call matvec(iiend-ii+1,kkend-kk+1,lda,
     $                           a(ii,kk) ,b(kk,j), c(ii,j) )
   30             continue
   40          continue
   50       continue
         else
            do 51 j = 1,p 
            call dgemv ( transa, m, n, alpha, a, lda, b( 1, j ), 1,
     $                   beta, c( 1, j ), 1 )
   51        continue
         end if
      else if( ( transa.eq.'n' ).and.( tb ) )then
         c( 1: m, 1: p ) = alpha*
     $                     matmul( a( 1: m, 1: n ),
     $                             transpose( b( 1: p, 1: n ) ) )     +
     $                     beta*c( 1: m, 1: p )
      else if( ta.and.tb )then
         c( 1: n, 1: p ) = alpha*
     $                     matmul( transpose( a( 1: m, 1: n ) ),
     $                             transpose( b( 1: p, 1: m ) ) )     +
     $                     beta*c( 1: n, 1: p )
      end if
*
      return
*
*     end of dgemm .
*
      end
End of dgemm.f
echo dgemv.f 1>&2
cat >dgemv.f <<'End of dgemv.f'
      SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
     $                   BETA, Y, INCY )
*     .. Scalar Arguments ..
      DOUBLE PRECISION   ALPHA, BETA
      INTEGER            INCX, INCY, LDA, M, N
      CHARACTER*1        TRANS
*     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
*     ..
*
*  Purpose
*  =======
*
*  DGEMV  performs one of the 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.
*
*  Parameters
*  ==========
*
*  TRANS  - CHARACTER*1.
*           On entry, TRANS specifies the operation to be performed as
*           follows:
*
*              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
*
*              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
*
*              TRANS = 'C' or 'c'   y := alpha*A'*x + beta*y.
*
*           Unchanged on exit.
*
*  M      - INTEGER.
*           On entry, M specifies the number of rows of the matrix A.
*           M must be at least zero.
*           Unchanged on exit.
*
*  N      - INTEGER.
*           On entry, N specifies the number of columns of the matrix A.
*           N must be at least zero.
*           Unchanged on exit.
*
*  ALPHA  - DOUBLE PRECISION.
*           On entry, ALPHA specifies the scalar alpha.
*           Unchanged on exit.
*
*  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
*           Before entry, the leading m by n part of the array A must
*           contain the matrix of coefficients.
*           Unchanged on exit.
*
*  LDA    - INTEGER.
*           On entry, LDA specifies the first dimension of A as declared
*           in the calling (sub) program. LDA must be at least
*           max( 1, m ).
*           Unchanged on exit.
*
*  X      - DOUBLE PRECISION array of DIMENSION at least
*           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
*           and at least
*           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
*           Before entry, the incremented array X must contain the
*           vector x.
*           Unchanged on exit.
*
*  INCX   - INTEGER.
*           On entry, INCX specifies the increment for the elements of
*           X. INCX must not be zero.
*           Unchanged on exit.
*
*  BETA   - DOUBLE PRECISION.
*           On entry, BETA specifies the scalar beta. When BETA is
*           supplied as zero then Y need not be set on input.
*           Unchanged on exit.
*
*  Y      - DOUBLE PRECISION array of DIMENSION at least
*           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
*           and at least
*           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
*           Before entry with BETA non-zero, the incremented array Y
*           must contain the vector y. On exit, Y is overwritten by the
*           updated vector y.
*
*  INCY   - INTEGER.
*           On entry, INCY specifies the increment for the elements of
*           Y. INCY must not be zero.
*           Unchanged on exit.
*
*
*  Level 2 Blas routine.
*
*  -- Written on 20-July-1986.
*     Jack Dongarra, Argonne National Lab.
*     Jeremy Du Croz, Nag Central Office.
*     Sven Hammarling, Nag Central Office.
*     Richard Hanson, Sandia National Labs.
*
*     This version for the Alliant FX series.
*     Written on 19-September-1986.
*     Sven Hammarling.
*
*     .. Parameters ..
      DOUBLE PRECISION   ONE         , ZERO
      PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
*     .. Local Scalars ..
      DOUBLE PRECISION   TEMP
      INTEGER            I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     .. External Subroutines ..
      EXTERNAL           XERBLA
*     .. Intrinsic Functions ..
      INTRINSIC          MAX
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      INFO = 0
      IF      ( .NOT.LSAME( TRANS, 'N' ).AND.
     $          .NOT.LSAME( TRANS, 'T' ).AND.
     $          .NOT.LSAME( TRANS, 'C' )      ) THEN
         INFO = 1
      ELSE IF ( M.LT.0 ) THEN
         INFO = 2
      ELSE IF ( N.LT.0 ) THEN
         INFO = 3
      ELSE IF ( LDA.LT.MAX(1,M) ) THEN
         INFO = 6
      ELSE IF ( INCX.EQ.0 ) THEN
         INFO = 8
      ELSE IF ( INCY.EQ.0 ) THEN
         INFO = 11
      END IF
      IF( INFO.NE.0 )THEN
         CALL XERBLA( 'DGEMV ', INFO )
         RETURN
      END IF
*
*     Quick return if possible.
*
      IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
     $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
     $   RETURN
*
*     Set LENX and LENY, the lengths of the vectors x and y.
*
      IF( LSAME( TRANS, 'N' ) )THEN
         LENX = N
         LENY = M
      ELSE
         LENX = M
         LENY = N
      END IF
*
*     Start the operations. In this version the elements of A are
*     accessed sequentially with one pass through A.
*
*     First form  y := beta*y  and set up the start points in X and Y if
*     the increments are not both unity.
*
      IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
         IF( BETA.NE.ONE )THEN
            IF( BETA.EQ.ZERO )THEN
               DO 10, I = 1, LENY
                  Y( I ) = ZERO
   10          CONTINUE
            ELSE
               DO 20, I = 1, LENY
                  Y( I ) = BETA*Y( I )
   20          CONTINUE
            END IF
         END IF
      ELSE
         IF( INCX.GT.0 )THEN
            KX = 1
         ELSE
            KX = 1 - ( LENX - 1 )*INCX
         END IF
         IF( INCY.GT.0 )THEN
            KY = 1
         ELSE
            KY = 1 - ( LENY - 1 )*INCY
         END IF
         IF( BETA.NE.ONE )THEN
            IY = KY
            IF( BETA.EQ.ZERO )THEN
               DO 30, I = 1, LENY
                  Y( IY ) = ZERO
                  IY      = IY   + INCY
   30          CONTINUE
            ELSE
               DO 40, I = 1, LENY
                  Y( IY ) = BETA*Y( IY )
                  IY      = IY           + INCY
   40          CONTINUE
            END IF
         END IF
      END IF
      IF( ALPHA.EQ.ZERO )
     $   RETURN
      IF( LSAME( TRANS, 'N' ) )THEN
*
*        Form  y := alpha*A*x + y.
*
         IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
            J = MOD( N, 2 )
            IF( J.EQ.1 )THEN
               Y( 1: M ) = Y( 1: M ) + ALPHA*X( 1 )*A( 1: M, 1 )
            END IF
            J = MOD( N, 4 )
            IF( ( J.EQ.2 ).OR.( J.EQ.3 ) )THEN
               Y( 1: M ) = Y( 1: M ) +
     $                     ALPHA*( X( J - 1 )*A( 1: M, J - 1 ) +
     $                             X( J     )*A( 1: M, J     )   )
            END IF
            J = MOD( N, 8 )
            IF( ( J.EQ.4 ).OR.( J.EQ.5 ).OR.
     $          ( J.EQ.6 ).OR.( J.EQ.7 ) )THEN
               Y( 1: M ) = Y( 1: M ) +
     $                     ALPHA*( X( J - 3 )*A( 1: M, J - 3 ) +
     $                             X( J - 2 )*A( 1: M, J - 2 ) +
     $                             X( J - 1 )*A( 1: M, J - 1 ) +
     $                             X( J     )*A( 1: M, J     )   )
            END IF
            JSTART = J + 8
            DO 60, J = JSTART, N, 8
C               IF( X( J ).NE.ZERO )THEN
C                  TEMP = ALPHA*X( J )
C                  DO 50, I = 1, M
C                     Y( I ) = Y( I ) + TEMP*A( I, J )
C   50             CONTINUE
                   Y( 1: M ) = Y( 1: M ) +
     $                         ALPHA*( X( J - 7 )*A( 1: M, J - 7 ) +
     $                                 X( J - 6 )*A( 1: M, J - 6 ) +
     $                                 X( J - 5 )*A( 1: M, J - 5 ) +
     $                                 X( J - 4 )*A( 1: M, J - 4 ) +
     $                                 X( J - 3 )*A( 1: M, J - 3 ) +
     $                                 X( J - 2 )*A( 1: M, J - 2 ) +
     $                                 X( J - 1 )*A( 1: M, J - 1 ) +
     $                                 X( J     )*A( 1: M, J     )   )
C               END IF
C            Y( 1: M ) = Y( 1: M ) +
C     $                  ALPHA*matmul( A( 1: M, 1: N ), X( 1: N ) )
   60       CONTINUE
         ELSE
            JX = KX
            DO 80, J = 1, N
               IF( X( JX ).NE.ZERO )THEN
                  TEMP = ALPHA*X( JX )
                  IY   = KY
                  DO 70, I = 1, M
                     Y( IY ) = Y( IY ) + TEMP*A( I, J )
                     IY      = IY      + INCY
   70             CONTINUE
               END IF
               JX = JX + INCX
   80       CONTINUE
         END IF
      ELSE
*
*        Form  y := alpha*A'*x + y.
*
         IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
            DO 100, J = 1, N
C               TEMP = ZERO
C               DO 90, I = 1, M
C                  TEMP = TEMP + A( I, J )*X( I )
C   90          CONTINUE
C               Y( J ) = Y( J ) + ALPHA*TEMP
               Y( J ) = Y( J ) + ALPHA*dotproduct( A( 1: M, J ),
     $                                             X( 1: M ) )
  100       CONTINUE
         ELSE
            JY = KY
            DO 120, J = 1, N
               TEMP = ZERO
               IX   = KX
               DO 110, I = 1, M
                  TEMP = TEMP + A( I, J )*X( IX )
                  IX   = IX   + INCX
  110          CONTINUE
               Y( JY ) = Y( JY ) + ALPHA*TEMP
               JY      = JY      + INCY
  120       CONTINUE
         END IF
      END IF
*
      RETURN
*
*     End of DGEMV .
*
      END
End of dgemv.f
echo lsame.f 1>&2
cat >lsame.f <<'End of lsame.f'
      LOGICAL FUNCTION LSAME ( CA, CB )
*     .. Scalar Arguments ..
      CHARACTER*1            CA, CB
*     ..
*
*  Purpose
*  =======
*
*  LSAME  tests if CA is the same letter as CB regardless of case.
*  CB is assumed to be an upper case letter. LSAME returns .TRUE. if
*  CA is either the same as CB or the equivalent lower case letter.
*
*  N.B. This version of the routine is only correct for ASCII code.
*       Installers must modify the routine for other character-codes.
*
*       For EBCDIC systems the constant IOFF must be changed to -64.
*       For CDC systems using 6-12 bit representations, the system-
*       specific code in comments must be activated.      
*
*  Parameters
*  ==========
*
*  CA     - CHARACTER*1
*  CB     - CHARACTER*1
*           On entry, CA and CB specify characters to be compared.
*           Unchanged on exit.
*
*
*  Auxiliary routine for Level 2 Blas.
*
*  -- Written on 20-July-1986
*     Richard Hanson, Sandia National Labs.
*     Jeremy Du Croz, Nag Central Office.
*
*     .. Parameters ..
      INTEGER                IOFF
      PARAMETER            ( IOFF=32 )
*     .. Intrinsic Functions ..
      INTRINSIC              ICHAR
*     .. Executable Statements ..
*
*     Test if the characters are equal
*
      LSAME = CA .EQ. CB
*
*     Now test for equivalence
*
      IF ( .NOT.LSAME ) THEN
         LSAME = ICHAR(CA) - IOFF .EQ. ICHAR(CB)
      END IF
*
      RETURN
*
*  The following comments contain code for CDC systems using 6-12 bit
*  representations.
*
*     .. Parameters ..
*     INTEGER                ICIRFX
*     PARAMETER            ( ICIRFX=62 )
*     .. Scalar Arguments ..
*     CHARACTER*1            CB
*     .. Array Arguments ..
*     CHARACTER*1            CA(*)
*     .. Local Scalars ..
*     INTEGER                IVAL
*     .. Intrinsic Functions ..
*     INTRINSIC              ICHAR, CHAR
*     .. Executable Statements ..
*
*     See if the first character in string CA equals string CB.
*
*     LSAME = CA(1) .EQ. CB .AND. CA(1) .NE. CHAR(ICIRFX)
*
*     IF (LSAME) RETURN
*
*     The characters are not identical. Now check them for equivalence.
*     Look for the 'escape' character, circumflex, followed by the
*     letter.
*
*     IVAL = ICHAR(CA(2))
*     IF (IVAL.GE.ICHAR('A') .AND. IVAL.LE.ICHAR('Z')) THEN
*        LSAME = CA(1) .EQ. CHAR(ICIRFX) .AND. CA(2) .EQ. CB
*     END IF
*
*     RETURN
*
*     End of LSAME.
*
      END
End of lsame.f
echo main.f 1>&2
cat >main.f <<'End of main.f'
      program            main
*
*     Jack Dongarra and Dan Sorensen, 1/26/87.
*     Argonne National Laboratory
*
      integer            m, n, p, maxr, maxc
      parameter        ( maxr = 513, maxc = 513 )
      double precision   a( maxr, maxc ), b( maxr, maxc )
      double precision   c( maxr, maxc ), cc( maxr, maxc )
      double precision   hrcdelta
      real               temp( 2 )
      integer            t3b( 2 ), t3e( 2 ), t3d( 2 )
      character*1        ta
*
c      write(6,*)' input ta,m,n,p'
c      read( 5, * )ta, m, n, p
      write(6,*)' input m,n,p'
      read( 5, * ) m, n, p
      ta = 'n'
*     repeat ...
c      m  = 512
   10 continue
c      n  = 512
c   20 continue
c      p  = 512
c   30 continue
         call matgen( ta, maxr, maxr, m, n, p, a, b )
         call zero(maxr,m,p,c)
         t1b = second( t )
         call dgemm( ta, 'n', m, n, p, 1.0d0, a, maxr, b, maxr,
     $               0.0d0, c, maxr )
         t1e = second( t )
         t1d = second( t )
         call zero(maxr,m,p,c)
         t2b = etime ( temp )
         call dgemm( ta, 'n', m, n, p, 1.0d0, a, maxr, b, maxr,
     $               0.0d0, c, maxr )
         t2e = etime ( temp )
         t2d = etime ( temp )
         call zero(maxr,m,p,c)
         call hrcget ( t3b )
         call dgemm( ta, 'n', m, n, p, 1.0d0, a, maxr, b, maxr,
     $               0.0d0, c, maxr )
         call hrcget ( t3e )
         call hrcget ( t3d )
         t1  = ( t1e - t1b ) - ( t1d - t1e )
         t2  = ( t2e - t2b ) - ( t2d - t2e )
         t3  = hrcdelta( t3b, t3e ) - hrcdelta( t3e, t3d )
         t   = ( ( 2.0*m )*n )*p
         if( t1.gt.0.0e0 )then
            mflop1 = t/( 1.0e6*t1 )
         else
            mflop1 = -1
         end if
         if( t2.gt.0.0e0 )then
            mflop2 = t/( 1.0e6*t2 )
         else
            mflop2 = -1
         end if
         if( t3.gt.0.0e0 )then
            mflop3 = t/( 1.0e6*t3 )
         else
            mflop3 = -1
         end if
         write( 6, 99 )m, n, p, mflop1, mflop2, mflop3
   99    format( ' m = ', i4, '  n = ', i4, '  p = ', i4/
     $           '                               mflop1 = ', i4/
     $           '                               mflop2 = ', i4/
     $           '                               mflop3 = ', i4/ )
         call matgen( ta, maxr, maxr, m, n, p, a, b )
         call check( ta.eq.'t',maxr, maxr, maxr, m, n, p, a, b, c, cc )
c         read( 5, * )ta, m, n, p
      write(6,*)' input m,n,p'
         read( 5, * ) m, n, p
c         p = p/2
c         if( p.gt.0 )go to 30
c         n = n/2
c         if( n.gt.0 )go to 20
c         m = m/2
c         if( m.gt.0 )go to 10
         if(        ( m.gt.0 ).and.( n.gt.0 ).and.( p.gt.0 ) )
     $      go to 10
*     ... until( ( m.gt.0 ).and.( n.gt.0 ).and.( p.gt.0 ) )
      stop
      end
End of main.f
echo matgen.f 1>&2
cat >matgen.f <<'End of matgen.f'
      subroutine matgen(ta,lda,ldb,n1,n2,n3,a,b)
c
*
*     Jack Dongarra and Dan Sorensen, 1/26/87.
*     Argonne National Laboratory
*
      character*1 ta
      integer lda,ldb,n1,n2,n3
      double precision a(lda,*),b(ldb,*)
c
      init = 3125
      do 20 j = 1,n2
         do 10 i = 1,n1
            init = mod(3125*init,65536)
            a(i,j) = (init - 32768.0)/16384.0
   10    continue
   20 continue
      if( ta.eq.'n' )then
         do 40 j = 1,n3
            do 30 i = 1,n2
               init = mod(3125*init,65536)
               b(i,j) = (init - 32768.0)/16384.0
   30       continue
   40    continue
      else
         do 60 j = 1,n3
            do 50 i = 1,n1
               init = mod(3125*init,65536)
               b(i,j) = (init - 32768.0)/16384.0
   50       continue
   60    continue
      end if
      return
      end
End of matgen.f
echo xerbla.f 1>&2
cat >xerbla.f <<'End of xerbla.f'
      SUBROUTINE XERBLA( SRNAME, INFO )
      INTEGER            INFO
      CHARACTER*6        SRNAME
*
*     This is a special version of XERBLA to be used only as part of
*     the test program for testing error exits from the Level 2 BLAS
*     routines.
*
*     XERBLA  is an error handler for the Level 2 BLAS routines.
*
*     It is called by the Level 2 BLAS routines if an input parameter is
*     invalid.
*
      LOGICAL OK, LERR
      CHARACTER*6 SRNAMT
      COMMON / INFOC / INFOT, NOUT, OK, LERR
      COMMON / SRNAMC / SRNAMT
      LERR = .TRUE.
      IF (INFO.NE.INFOT) THEN
         WRITE (NOUT,99998) INFO, INFOT
         OK = .FALSE.
      ENDIF
      IF (SRNAME.NE.SRNAMT) THEN
         WRITE (NOUT,99997) SRNAME, SRNAMT
         OK = .FALSE.
      END IF
      RETURN
*
99998 FORMAT (' ******* XERBLA WAS CALLED WITH INFO = ', I6,
     $        ' INSTEAD OF ',I2,' *******')
99997 FORMAT (' ******* XERBLA WAS CALLED WITH SRNAME = ', A6,
     $        ' INSTEAD OF ',A6,' *******')
      END
End of xerbla.f
echo zero.f 1>&2
cat >zero.f <<'End of zero.f'
      subroutine zero(ldc,n1,n3,c)
*
*     Jack Dongarra and Dan Sorensen, 1/26/87.
*     Argonne National Laboratory
*
      integer ldc,n1,n3
      double precision c(ldc,*)
      do 20 j = 1,n3
         do 10 i = 1,n1
            c(i,j) = 0.0
   10    continue
   20 continue
      return
      end
End of zero.f
echo mvmul.s 1>&2
cat >mvmul.s <<'End of mvmul.s'


|                        Machine Code Listing of matvec.f 
|
|                       0 4 8   12 16 20
|     subroutine matvec(m,n,lda,a ,x, y )
|     double precision a(lda,1),x(lda,1),y(lda,1)
|     do 100 j = 1,n
|           y(1:m,1:4) = y(1:m,1:4) + a(1:m,j)*x(j,1:4)
| 100 continue
|     return
|     end
|     
|
|     Jack Dongarra and Dan Sorensen, 1/26/87.
|     Argonne National Laboratory
|
       .data
       .bss
       .text
_BTEXT:

| Literal Pool

       
       .globl _matvec_
_matvec_:                         |> 0000 
       link a6,#-40               |  grab stack space
                                  | LINE 2
       movl a0@(8),a2             |  put address of lda in a2
                                  | LINE 5
       movl a0@(4),a1             |  put address of n in a1
                                  | LINE 2
       movl a0@(20),a5            |  put address of y in a5
       movl a0@(12),a4            |  put address of a in a4
       movl a2@,d2                |  put value of lda in d2
                                  | LINE 5
       movl a0@,a3                |  put value of m in a3
       movl a3@,d4                |  into d4 to set vec length
       movl a0@(16),a3            |  move address of x into a3
       moveq #0,d0                |  d0 will point to row blocks of a     
       moveq #-1,d6               |  set vector mask to all ones
       moveq #1,d5                |  set stride to one
matvec.label_L1:
       movl a1@,d7                |  put value of n in d7
       subql #1,d7                |  initialize loop counter
       moveq #1,d3                |  d3 indexes the x vector
       movl  d0,d1                |  d1 indexes the columns of a
       vmoved a5@(0:b)[d0:l:d],v1 |  move y1 into v1     
       addl d2,d1
       vmoved a5@(0:b)[d1:l:d],v2 |  move y2 into v2     
       addl d2,d1
       vmoved a5@(0:b)[d1:l:d],v3 |  move y3 into v3     
       addl d2,d1
       vmoved a5@(0:b)[d1:l:d],v4 |  move y4 into v4     
       movl  d0,d1                |  d1 indexes the columns of a
matvec.label_LE:
                                  | LINE 6
       fmoved a3@(-8:b)[d3:l:d],fp0
                                  | get jth component xj of x into fp0
       vmuadd fp0,a4@(0:b)[d1:l:d],v1,v1                
                                  | mult col of a by xj and add to result
       addl d2,d3
       fmoved a3@(-8:b)[d3:l:d],fp0
       vmuadd fp0,a4@(0:b)[d1:l:d],v2,v2                
                                  | mult col of a by xj and add to result
       addl d2,d3
       fmoved a3@(-8:b)[d3:l:d],fp0
       vmuadd fp0,a4@(0:b)[d1:l:d],v3,v3                
                                  | mult col of a by xj and add to result
       addl d2,d3
       fmoved a3@(-8:b)[d3:l:d],fp0
       vmuadd fp0,a4@(0:b)[d1:l:d],v4,v4                
                                  | mult col of a by xj and add to result
       subl d2,d3
       subl d2,d3
       subl d2,d3
                                  | LINE 7
       addl d5,d3                 | increment x index
                                  | LINE 6
       addl d2,d1                 | increment column index
       subql #1,d7                | decrement loop counter
                                  | LINE 9
       bges matvec.label_LE       |> 006A  6CCC
       vmoved v1,a5@(0:b)[d0:l:d] | write result to y1
       movl  d0,d1                |  d1 indexes the columns of a
       addl d2,d1
       vmoved v2,a5@(0:b)[d1:l:d] | write result to y2
       addl d2,d1
       vmoved v3,a5@(0:b)[d1:l:d] | write result to y3
       addl d2,d1
       vmoved v4,a5@(0:b)[d1:l:d] | write result to y4
       moveq  #32,d7
       addl   d7,d0
       vcnt32 matvec.label_L1
       unlk a6                    |> 006C  4E5E
       rts                        |> 006E  4E75
End of mvmul.s
echo makefile 1>&2
cat >makefile <<'End of makefile'
FILES3 = check.o main.o matgen.o dgemm.o  dgemv.o lsame.o xerbla.o mvmul.o zero.o

run : $(FILES3)
	fortran -recursive -O -AS $(FILES3) /afs2/dongarra/lib.a -o run

.f.o : ; fortran -recursive -O -AS -c $*.f
End of makefile
