!*********************************************************
!*********************************************************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
