PC电路运算速度性能测试基准程序

!*********************************************************
!*********************************************************72
Program main

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

!c gprof_test, illustrate the use of the GPROF program

!                                    performance monitor

!********************************************************
  Double Precision a(1001, 1000), b(1000), x(1000)

  Double Precision time(6), cray, ops, total, norma, normx

  Double Precision resid, residn, eps, epslon

  Integer ipvt(1000)

  lda = 1001

!

!     this program was updated on 10/12/92 to correct a

!     problem with the random number generator. The previous

!     random number generator had a short period and produced

!     singular matrices occasionally.


!
  Call timestamp()


  Write (*, '(a)') ' '

  n = 1000

  cray = .056

  Write (6, 1)

  1 format('$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'//&
  '新年好！ 兔年吉祥!'/&
                '玉兔报喜，福星献瑞！'/&
                     '虎挟雄风辞旧岁，玉兔呈祥月中来!'//&
      '$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'/&
      '注意防疫，加油！'/&
        'Protect yourself well from the COVID19-Pandemics!'/&
   '$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'//&
             'Merry Christmas & Happy New Year!'//&
             'Wish you good health & prosperous!'//&
             'SINOMACH'//&
             'PotsyYZ, 周勇, GoldbergJr.JohnBkt'/&
             '广州    20221214'/&
             '$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'//&
' Please send the results of this run to:'//&


  ' Jack J. Dongarra'/&


  ' Computer Science Department'/&


  ' University of Tennessee'/&


  ' Knoxville, Tennessee 37996-1300'//&


  ' Fax: 615-974-8296'//&


  ' Internet: dongarra@cs.utk.edu' /)

  ops = (2.0D0*dfloat(n)**3)/3.0D0 + 2.0D0*dfloat(n)**2

!

  Call matgen(a, lda, n, b, norma)

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

!        you should replace the call to dgefa and dgesl

!        by calls to your linear equation solver.

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

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

!

  t1 = second()


  Call dgefa(a, lda, n, ipvt, info)

  time(1) = second() - t1

  t1 = second()

  Call dgesl(a, lda, n, ipvt, b, 0)

  time(2) = second() - t1

  total = time(1) + time(2)

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

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

!

!     compute a residual to verify results.

!

  Do i = 1, n


    x(i) = b(i)


  End Do



  Call matgen(a, lda, n, b, norma)

  Do i = 1, n

    b(i) = -b(i)

  End Do


  Call dmxpy(n, b, n, lda, x, a)

  resid = 0.0

  normx = 0.0

  Do i = 1, n

    resid = dmax1(resid, dabs(b(i)))

    normx = dmax1(normx, dabs(x(i)))

  End Do



  eps = epslon(1.0D0)

  residn = resid/(n*norma*normx*eps)

  Write (6, 40)

  Write (6, 50) residn, resid, eps, x(1), x(n)


!

  Write (6, 60) n

  Write (6, 70)

  70 format(6 x, 'factor', 5 x, 'solve', 6 x, 'total', 5 x, 'mflops', 7 x,&

  'unit',6x,'ratio')

!

  time(3)=total

  time(4)=ops/(1.0D6*total)

  time(5)=2.0D0/time(4)

  time(6)=total/cray

  Write(6,80) lda



  Write(6,110)(time(i),i=1,6)


  Write(6,*) ' end of tests -- this version dated 10/12/92'

!

  Write(*,'(a)') ' '


  Call timestamp()

  Stop

  40 Format('     norm. resid      resid           machep',&


  '         x(1)          x(n)')

  50 Format(1P5E16.8)

  60 Format(//'    times are reported for matrices of order ',I5)


  80 Format(' times for array with leading dimension of',I4)


  110 Format(6(1PE11.3))

End Program main




Subroutine matgen(a,lda,n,b,norma)

  Integer lda, n, init(4), i, j

  Double Precision a(lda,1), b(1), norma, random_value

!
  init(1)=1

  init(2)=2

  init(3)=3

  init(4)=1325

  norma=0.0

  Do j=1, n

    Do i=1, n

      a(i,j)=random_value(init)-.5

      norma=dmax1(dabs(a(i,j)),norma)

    End Do

  End Do


  Do i=1, n


    b(i)=0.0


  End Do


  Do j=1, n


    Do i=1, n


      b(i)=b(i)+a(i,j)


    End Do


  End Do


  Return


End Subroutine matgen




Subroutine dgefa(a,lda,n,ipvt,info)

  Integer lda, n, ipvt(1), info


  Double Precision a(lda,1)


!

!     dgefa factors a double precision matrix by gaussian elimination.

!

!     dgefa is usually called by dgeco, but it can be called

!     directly with a saving in time if  rcond  is not needed.

!     (time for dgeco) = (1 + 9/n)*(time for dgefa) .

!

!     on entry

!

!        a       double precision(lda, n)

!                the matrix to be factored.

!

!        lda     integer

!                the leading dimension of the array  a .

!

!        n       integer

!                the order of the matrix  a .

!


!     on return

!

!        a       an upper triangular matrix and the multipliers

!                which were used to obtain it.

!                the factorization can be written  a = l*u  where

!                l  is a product of permutation and unit lower

!                triangular matrices and  u  is upper triangular.

!

!        ipvt    integer(n)

!                an integer vector of pivot indices.

!

!        info    integer

!                = 0  normal value.

!                = k  if  u(k,k) .eq. 0.0 .  this is not an error

!                     condition for this subroutine, but it does

!                     indicate that dgesl or dgedi will divide by zero

!                     if called.  use  rcond  in dgeco for a reliable

!                     indication of singularity.

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

!     linpack. this version dated 08/14/78 .

!     cleve moler, university of new mexico, argonne national lab.

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



!
!     subroutines and functions

!     blas daxpy,dscal,idamax

!     internal variables

!

  Double Precision t

  Integer idamax, j, k, kp1, l, nm1

!

!     gaussian elimination with partial pivoting

!
  info=0


  nm1=n-1


  If(nm1<1) Goto 70


  Do k=1, nm1


    kp1=k+1

!

!        find l = pivot index

!

    l=idamax(n-k+1,a(k,k),1)+k-1


    ipvt(k)=l


!

!        zero pivot implies this column already triangularized

!

    If(a(l,k)==0.0D0) Goto 40

!

!           interchange if necessary

!

    If(l==k) Goto 10

    t=a(l,k)

    a(l,k)=a(k,k)

    a(k,k)=t


    10 Continue

!

!           compute multipliers

!

    t=-1.0D0/a(k,k)



    Call dscal(n-k,t,a(k+1,k),1)


!

!           row elimination with column indexing

!

    Do j=kp1, n


       t=a(l,j)


      If(l==k) Goto 20


      a(l,j)=a(k,j)


      a(k,j)=t


      20 Continue




      Call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)

    End Do


    Goto 50

    40 Continue


    info=k

    50 Continue

  End Do


  70 Continue


  ipvt(n)=n


  If(a(n,n)==0.0D0) info=n


  Return

End Subroutine dgefa





Subroutine dgesl(a,lda,n,ipvt,b,job)

  Integer lda, n, ipvt(1), job

  Double Precision a(lda,1), b(1)

!

!     dgesl solves the double precision system

!     a * x = b  or  trans(a) * x = b


!     using the factors computed by dgeco or dgefa.

!

!     on entry

!

!        a       double precision(lda, n)

!                the output from dgeco or dgefa.

!

!        lda     integer

!                the leading dimension of the array  a .

!

!        n       integer

!                the order of the matrix  a .

!

!        ipvt    integer(n)

!                the pivot vector from dgeco or dgefa.

!

!        b       double precision(n)

!                the right hand side vector.
!

!        job     integer

!                = 0         to solve  a*x = b ,

!                = nonzero   to solve  trans(a)*x = b  where

!                            trans(a)  is the transpose.

!

!     on return

!

!        b       the solution vector  x .

!

!     error condition

!

!        a division by zero will occur if the input factor contains a

!        zero on the diagonal.  technically this indicates singularity

!        but it is often caused by improper arguments or improper

!        setting of lda .  it will not occur if the subroutines are

!        called correctly and if dgeco has set rcond .gt. 0.0

!        or dgefa has set info .eq. 0 .

!

!     to compute  inverse(a) * c  where  c  is a matrix

!     with  p  columns

!           call dgeco(a,lda,n,ipvt,rcond,z)

!           if (rcond is too small) go to ...

!           do 10 j = 1, p

!              call dgesl(a,lda,n,ipvt,c(1,j),0)

!        10 continue

!


!     linpack. this version dated 08/14/78 .

!     cleve moler, university of new mexico, argonne national lab.

!

!     subroutines and functions

!

!     blas daxpy,ddot

!

!     internal variables

!
  Double Precision ddot, t

  Integer k, kb, l, nm1

!
  nm1=n-1

  If(job/=0) Goto 50

!

!        job = 0 , solve  a * x = b

!        first solve  l*y = b

!

  If(nm1<1) Goto 30


  Do k=1, nm1


    l=ipvt(k)


    t=b(l)


    If(l==k) Goto 10


    b(l)=b(k)


    b(k)=t


    10 Continue


    Call daxpy(n-k,t,a(k+1,k),1,b(k+1),1)


  End Do


  30 Continue


!

!        now solve  u*x = y

!

  Do kb=1, n

    k=n+1-kb

    b(k)=b(k)/a(k,k)

    t=-b(k)



    Call daxpy(k-1,t,a(1,k),1,b(1),1)


  End Do


  Goto 100


  50 Continue


!

!        job = nonzero, solve  trans(a) * x = b

!        first solve  trans(u)*y = b

!

  Do k=1, n

    t=ddot(k-1,a(1,k),1,b(1),1)

    b(k)=(b(k)-t)/a(k,k)

  End Do

!

!        now solve trans(l)*x = y

!

  If(nm1<1) Goto 90


  Do kb=1, nm1


    k=n-kb


    b(k)=b(k)+ddot(n-k,a(k+1,k),1,b(k+1),1)


    l=ipvt(k)


    If(l==k) Goto 70


    t=b(l)


    b(l)=b(k)


    b(k)=t


    70 Continue


  End Do



  90 Continue



  100 Continue


  Return

End Subroutine dgesl






Subroutine daxpy(n,da,dx,incx,dy,incy)

!

!     constant times a vector plus a vector.

!     uses unrolled loops for increments equal to one.


!     jack dongarra, linpack, 3/11/78.

!

  Double Precision dx(1), dy(1), da

  Integer i, incx, incy, ix, iy, m, mp1, n

!

  If(n<=0) Return

  If(da==0.0D0) Return

  If(incx==1 .And. incy==1) Goto 20

!

!        code for unequal increments or equal increments

!          not equal to 1

!

  ix=1


  iy=1


  If(incx<0) ix=(-n+1)*incx+1


  If(incy<0) iy=(-n+1)*incy+1


  Do i=1, n


    dy(iy)=dy(iy)+da*dx(ix)


    ix=ix+incx


    iy=iy+incy


  End Do


  Return


!

!        code for both increments equal to 1

!

!

!        clean-up loop

!

  20 m=mod(n,4)


  If(m==0) Goto 40


  Do i=1, m


    dy(i)=dy(i)+da*dx(i)


  End Do



  If(n<4) Return


  40 mp1=m+1


  Do i=mp1, n, 4


    dy(i)=dy(i)+da*dx(i)


    dy(i+1)=dy(i+1)+da*dx(i+1)


    dy(i+2)=dy(i+2)+da*dx(i+2)


    dy(i+3)=dy(i+3)+da*dx(i+3)


  End Do


  Return


End Subroutine daxpy



Double Precision Function ddot(n,dx,incx,dy,incy)


!

!     forms the dot product of two vectors.

!     uses unrolled loops for increments equal to one.

!     jack dongarra, linpack, 3/11/78.

!

  Double Precision dx(1), dy(1), dtemp

  Integer i, incx, incy, ix, iy, m, mp1, n

!

  ddot=0.0D0


  dtemp=0.0D0


  If(n<=0) Return


  If(incx==1 .And. incy==1) Goto 20


!

!        code for unequal increments or equal increments

!          not equal to 1

!

  ix=1

  iy=1

  If(incx<0) ix=(-n+1)*incx+1

  If(incy<0) iy=(-n+1)*incy+1


  Do i=1, n



    dtemp=dtemp+dx(ix)*dy(iy)



    ix=ix+incx



    iy=iy+incy


  End Do



  ddot=dtemp


  Return


!

!        code for both increments equal to 1


!        clean-up loop

!


  20 m=mod(n,5)


  If(m==0) Goto 40


  Do i=1, m


    dtemp=dtemp+dx(i)*dy(i)

  End Do




  If(n<5) Goto 60



  40 mp1=m+1



  Do i=mp1, n, 5



    dtemp=dtemp+dx(i)*dy(i)+dx(i+1)*dy(i+1)+dx(i+2)*dy(i+2)+dx(i+3)*dy(i+3)&



    +dx(i+4)*dy(i+4)



  End Do



  60 ddot=dtemp


  Return


End Function ddot





Subroutine dscal(n,da,dx,incx)

!

!     scales a vector by a constant.

!     uses unrolled loops for increment equal to one.

!     jack dongarra, linpack, 3/11/78.

!


  Double Precision da, dx(1)


  Integer i, incx, m, mp1, n, nincx

!


  If(n<=0) Return




  If(incx==1) Goto 20


!

!        code for increment not equal to 1


!

  nincx=n*incx



  Do i=1, nincx, incx



    dx(i)=da*dx(i)



  End Do



  Return


!

!        code for increment equal to 1



!        clean-up loop


!

  20 m=mod(n,5)



  If(m==0) Goto 40



  Do i=1, m



    dx(i)=da*dx(i)



  End Do



  If(n<5) Return



  40 mp1=m+1



  Do i=mp1, n, 5



    dx(i)=da*dx(i)



    dx(i+1)=da*dx(i+1)



    dx(i+2)=da*dx(i+2)



    dx(i+3)=da*dx(i+3)



    dx(i+4)=da*dx(i+4)



  End Do



  Return



End Subroutine dscal




Integer Function idamax(n,dx,incx)

!

!     finds the index of element having max. dabsolute value.

!     jack dongarra, linpack, 3/11/78.

!

  Double Precision dx(1), dmax


  Integer i, incx, ix, n

!

  idamax=0

  If(n<1) Return


  idamax=1


  If(n==1) Return


  If(incx==1) Goto 20

!

!        code for increment not equal to 1

!

  ix=1


  dmax=dabs(dx(1))


  ix=ix+incx


  Do i=2, n


    If(dabs(dx(ix))<=dmax) Goto 5


    idamax=i


    dmax=dabs(dx(ix))


    5 ix=ix+incx


  End Do


  Return


!

!        code for increment equal to 1

!

  20 dmax=dabs(dx(1))


  Do i=2, n


    If(dabs(dx(i))<=dmax) Goto 30


    idamax=i


    dmax=dabs(dx(i))


  30 End Do


  Return

End Function idamax





Double Precision Function epslon(x)



  Double Precision x

!

!     estimate unit roundoff in quantities of size x.

!

  Double Precision a, b, c, eps

!

!     this program should function properly on all systems

!     satisfying the following two assumptions,

!        1.  the base used in representing dfloating point

!            numbers is not a power of three.

!        2.  the quantity  a  in statement 10 is represented to

!            the accuracy used in dfloating point variables

!            that are stored in memory.

!     the statement number 10 and the go to 10 are intended to


!     force optimizing compilers to generate code satisfying


!     assumption 2.


!     under these assumptions, it should be true that,


!            a  is not exactly equal to four-thirds,


!            b  has a zero for its last bit or digit,


!            c  is not exactly equal to one,


!            eps  measures the separation of 1.0 from


!                 the next larger dfloating point number.


!     the developers of eispack would appreciate being informed


!     about any systems where these assumptions do not hold.

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

!     this routine is one of the auxiliary routines used by eispack iii

!     to avoid machine dependencies.

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

!     this version dated 4/6/83.

!

  a=4.0D0/3.0D0


  10 b=a-1.0D0


  c=b+b+b


  eps=dabs(c-1.0D0)


  If(eps==0.0D0) Goto 10


  epslon=eps*dabs(x)


  Return


End Function epslon





Subroutine mm(a,lda,n1,n3,b,ldb,n2,c,ldc)


  Double Precision a(lda,*), b(ldb,*), c(ldc,*)

!

!   purpose:

!     multiply matrix b times matrix c and store the result in matrix a.

!

!   parameters:

!

!     a double precision(lda,n3), matrix of n1 rows and n3 columns

!

!     lda integer, leading dimension of array a

!

!     n1 integer, number of rows in matrices a and b

!

!     n3 integer, number of columns in matrices a and c

!

!     b double precision(ldb,n2), matrix of n1 rows and n2 columns

!

!     ldb integer, leading dimension of array b

!

!     n2 integer, number of columns in matrix b, and number of rows in


!         matrix c


!


!     c double precision(ldc,n3), matrix of n2 rows and n3 columns

!

!     ldc integer, leading dimension of array c

!

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

!
  Do j=1, n3



    Do i=1, n1



      a(i,j)=0.0



    End Do



    Call dmxpy(n2,a(1,j),n1,ldb,c(1,j),b)


  End Do


!

  Return

End Subroutine mm




Subroutine dmxpy(n1,y,n2,ldm,x,m)


  Double Precision y(*), x(*), m(ldm,*)


!

!   purpose:


!     multiply matrix m times vector x and add the result to vector y.


!


!   parameters:


!


!     n1 integer, number of elements in vector y, and number of rows in


!         matrix m


!


!     y double precision(n1), vector of length n1 to which is added


!         the product m*x


!


!     n2 integer, number of elements in vector x, and number of columns


!         in matrix m


!


!     ldm integer, leading dimension of array m


!


!     x double precision(n2), vector of length n2


!


!     m double precision(ldm,n2), matrix of n1 rows and n2 columns


!

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

!

!   cleanup odd vector

!

  j=mod(n2,2)


  If(j>=1) Then


    Do i=1, n1


      y(i)=(y(i))+x(j)*m(i,j)


    End Do


  End If


!


!   cleanup odd group of two vectors


!


  j=mod(n2,4)



  If(j>=2) Then



    Do i=1, n1


      y(i)=((y(i))+x(j-1)*m(i,j-1))+x(j)*m(i,j)


    End Do


  End If


!

!   cleanup odd group of four vectors


!


  j=mod(n2,8)


  If(j>=4) Then


    Do i=1, n1



      y(i)=((((y(i))+x(j-3)*m(i,j-3))+x(j-2)*m(i,j-2))+x(j-1)*m(i,j-1))+&



      x(j)*m(i,j)


    End Do



  End If


!


!   cleanup odd group of eight vectors


!


  j=mod(n2,16)


  If(j>=8) Then


    Do i=1, n1


      y(i)=((((((((y(i))+x(j-7)*m(i,j-7))+x(j-6)*m(i,j-6))+x(j-5)*m(i,j-5))+&


      x(j-4)*m(i,j-4))+x(j-3)*m(i,j-3))+x(j-2)*m(i,j-2))+x(j-1)*m(i,j-1))+&


      x(j)*m(i,j)


    End Do

  End If

!

!   main loop - groups of sixteen vectors

!

  jmin=j+16


  Do j=jmin, n2, 16


    Do i=1, n1


      y(i)=((((((((((((((((y(i))&


      +x(j-15)*m(i,j-15))+x(j-14)*m(i,j-14))&


      +x(j-13)*m(i,j-13))+x(j-12)*m(i,j-12))&




      +x(j-11)*m(i,j-11))+x(j-10)*m(i,j-10))&


      +x(j-9)*m(i,j-9))+x(j-8)*m(i,j-8))&


      +x(j-7)*m(i,j-7))+x(j-6)*m(i,j-6))&


      +x(j-5)*m(i,j-5))+x(j-4)*m(i,j-4))&


      +x(j-3)*m(i,j-3))+x(j-2)*m(i,j-2))&


      +x(j-1)*m(i,j-1))+x(j)*m(i,j)


    End Do

  End Do


  Return

End Subroutine dmxpy




Double Precision Function random_value(iseed)

!

!     modified from the LAPACK auxiliary routine 10/12/92 JD

!  -- LAPACK auxiliary routine (version 1.0) --

!     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,

!     Courant Institute, Argonne National Lab, and Rice University

!     February 29, 1992

!

!     .. Array Arguments ..

  Integer iseed(4)

!     ..

!

!  Purpose

!  =======

!

!  DLARAN returns a random real number from a uniform (0,1)


!  distribution.

!

!  Arguments


!  =========

!

!  ISEED   (input/output) INTEGER array, dimension (4)


!          On entry, the seed of the random number generator; the array


!          elements must be between 0 and 4095, and ISEED(4) must be


!          odd.


!          On exit, the seed is updated.


!


!  Further Details


!  ===============


!


!  This routine uses a multiplicative congruential method with modulus


!  2**48 and multiplier 33952834046453 (see G.S.Fishman,



!  'Multiplicative congruential random number generators with modulus



!  2**b: an exhaustive analysis for b = 32 and a partial analysis for



!  b = 48', Math. Comp. 189, pp 331-344, 1990).


!


!  48-bit integers are stored in 4 integer array elements with 12 bits


!  per element. Hence the routine is portable across machines with


!  integers of 32 bits or more.


!


!     .. Parameters ..


  Integer m1, m2, m3, m4


  Parameter(m1=494,m2=322,m3=2508,m4=2549)


  Double Precision one


  Parameter(one=1.0D+0)



  Integer ipw2



  Double Precision r



  Parameter(ipw2=4096,r=one/ipw2)


!     ..


!     .. Local Scalars ..


  Integer it1, it2, it3, it4


!     ..


!     .. Intrinsic Functions ..


  Intrinsic dble, mod


!     ..


!     .. Executable Statements ..


!


!     multiply the seed by the multiplier modulo 2**48


!

  it4=iseed(4)*m4


  it3=it4/ipw2


  it4=it4-ipw2*it3


  it3=it3+iseed(3)*m4+iseed(4)*m3


  it2=it3/ipw2


  it3=it3-ipw2*it2


  it2=it2+iseed(2)*m4+iseed(3)*m3+iseed(4)*m2


  it1=it2/ipw2


  it2=it2-ipw2*it1


  it1=it1+iseed(1)*m4+iseed(2)*m3+iseed(3)*m2+iseed(4)*m1



  it1=mod(it1,ipw2)


!

!     return updated seed


!

  iseed(1)=it1



  iseed(2)=it2



  iseed(3)=it3



  iseed(4)=it4


!

!     convert 48-bit integer to a real number in the interval (0,1)

!

  random_value=r*(dble(it1)+r*(dble(it2)+r*(dble(it3)+r*(dble(it4)))))

  Return

!

!     End of RAN

!

End Function random_value




Subroutine timestamp()

!*********************************************************************72

!

!c TIMESTAMP prints out the current YMDHMS date as a timestamp.

!

!  Parameters:

!

!    None

!
  Implicit None

  Character*(8) ampm

  Integer d

  Character*(8) date

  Integer h

  Integer m

  Integer mm


  Character*(9) month(12)


  Integer n


  Integer s


  Character*(10) time


  Integer y


  Save month


  Data month/'January  ', 'February ', 'March    ', 'April    ', 'May      ', &


  'June     ', 'July     ', 'August   ', 'September', 'October  ', &


  'November ', 'December '/



  Call date_and_time(date,time)


  Read(date,'(i4,i2,i2)') y, m, d


  Read(time,'(i2,i2,i2,1x,i3)') h, n, s, mm



  If(h<12) Then



    ampm='AM'


  Else If(h==12) Then


    If(n==0 .And. s==0) Then


      ampm='Noon'



    Else



      ampm='PM'



    End If



  Else



    h=h-12



    If(h<12) Then



      ampm='PM'



    Else If(h==12) Then



      If(n==0 .And. s==0) Then



        ampm='Midnight'



      Else



        ampm='AM'



      End If



    End If



  End If


  Write(*,'(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)') d,&


  month(m), y, h, ':', n, ':', s, '.', mm, ampm


  Return


End Subroutine timestamp