完善资料让更多小伙伴认识你,还能领取20积分哦, 立即完善>
!********************************************************* !*********************************************************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 |
|
相关推荐 |
|
只有小组成员才能发言,加入小组>>
3448 浏览 4 评论
3487 浏览 1 评论
4348 浏览 1 评论
4940 浏览 8 评论
功率很小的反激电源理论上低压输入时采用CCM模式原边电流的有效值会小,但是看公式好像CCM模式更大一些
5254 浏览 1 评论
1202浏览 0评论
小黑屋| 手机版| Archiver| 电子发烧友 ( 湘ICP备2023018690号 )
GMT+8, 2024-11-21 16:47 , Processed in 1.040639 second(s), Total 50, Slave 37 queries .
Powered by 电子发烧友网
© 2015 bbs.elecfans.com
关注我们的微信
下载发烧友APP
电子发烧友观察
版权所有 © 湖南华秋数字科技有限公司
电子发烧友 (电路图) 湘公网安备 43011202000918 号 电信与信息服务业务经营许可证:合字B2-20210191 工商网监 湘ICP备2023018690号