FORTRAN 90 and Python KRISS http://ihlee.kriss.re.kr/compphys/f90module.htm http://ihlee.kriss.re.kr/compphys/pfphysics.htm Stay hungry, stay foolish.
Language is the house of the truth of Being. -Martin Heidegger.
Computer programming (often shortened to programming or coding) is the process of writing, testing, debugging/troubleshooting, and maintaining the source code of computer programs. A programming language is an artificial language that can be used to control the behavior of a machine, particularly a computer. A prominent purpose of programming languages is to provide instructions to a computer. Efficiency, reliability, robustness, portability, readability :!,. Fortran 90 and Python. -Martin Heidegger 1.. 2.. 3.. debugging Science..
Debugging is a very important task for every programmer, because an erroneous program is often useless. Bugs - Compile-time Errors Bugs - Run-time Errors 1945 9 9 bug; `debugging`..
Why FORTRAN 90? FORTRAN 90 1.. ( ;.) 2.. (, 113, ) 3.. (,, ) 4.. (,, ) 5.. ( 77.) ()?: (object), ( /...)... 90.... 1. data abstraction -- user-defined types 2. data hiding -- private and public attributes 3. encapsulation -- modules and data hiding facilities 4. inheritance and extensibility -- generic procedures 5. polymorphism -- generic overloading 6. reusability -- modules, ( ) John Backus
!234567890 program HW implicit none write(6,*) 'Hello World!' 77. C c *! Comment F90.! comment. end program hw ihlee@cetus0:~/lecture_programming$ a.out Hello World! ihlee@cetus0:~/lecture_programming$ 90 : http://incredible.egloos.com/3044862
. FORTRNA 90.
abc=1.0d0 abc 1.d0. 1.e0 (single precision) vs 1.d0 (double precision) 1.d0. (,,,,,,.) () Implicit Typing Undeclared variables have an implicit type, if first letter is I, J, K, L, M or N then type is INTEGER; any other letter then type is REAL s. Operator precedence pi*radius**2 = pi*(radius**2) x*y/z+a= (x*y/z)+a
, Intrinsic Logical Operations A LOGICAL or boolean expression returns a.true. /.FALSE. result. The following are valid with LOGICAL operands,.not. --.TRUE. if operand is.false...and. --.TRUE. if both operands are.true.;.or. --.TRUE. if at least one operand is.true.;.eqv. --.TRUE. if both operands are the same;.neqv. --.TRUE. if both operands are different. For example, if T is.true. and F is.false. FORTRAN 77 do while.. while (logical expr) do statements enddo do while (logical expr) statements enddo integer n n = 1 10 if (n.le. 100) then n = 2*n write (*,*) n goto 10 endif
IF (I > 17) THEN Print*, "I > 17" ELSE Print*, "I <= 17" END IF IF (I > 17) THEN Print*, "I > 17" ELSEIF (I == 17) THEN Print*, "I == 17" ELSE Print*, "I < 17" END IF outa: IF (a.ne. 0) THEN PRINT*, "a /= 0" IF (c.ne. 0) THEN PRINT*, "a /= 0 AND c /= 0" ELSE PRINT*, "a /= 0 BUT c == 0" ENDIF ELSEIF (a.gt. 0) THEN outa PRINT*, "a > 0" ELSE outa PRINT*, "a must be < 0" ENDIF outa [name:] IF (expression) THEN block ENDIF [name] [name:] IF (expression) THEN block1 ELSE [name] block2 ENDIF [name] elseif.. IF ( n == 1 ) THEN discount = 0.0d0 ELSEIF ( n <=5 ) THEN discount = 0.1d0 ELSEIF ( n>5.and. n <=10) THEN discount = 0.15d0 ELSE discount = 0.25d0 ENDIF
INTEGER :: month season: SELECT CASE( month ) CASE(4,5) WRITE(*,*) `Spring' CASE(6,7) WRITE(*,*) `Summer' CASE(8:10) WRITE(*,*) `Autumn' CASE(11,1:3,12) WRITE(*,*) `Winter' CASE DEFAULT WRITE(*,*) `not a month' END SELCET season Conditional statements IF statement IF construct CASE construct Repetition DO constructs GOTO outer: DO... inner: DO k=1,5 IF ( a(k)==0 ) CYCLE IF ( a(k)<0 ) EXIT outer... END DO inner END DO outer Exit from loop with EXIT. Begin another iteration with CYCLE. By default EXIT and CYCLE apply to the inner-most loop containing the statement, but can refer to a specific, named loop. INTEGER :: x... IF ( x.lt.1 ) GOTO 10... 10 x = 1
, integer i,j,k,nn integer mm real*8 abc real*8 cde(100),vec(3),xaa(3,3) real*8, allocatable :: dd(:,:), ee(:) integer, allocatable :: kkx(:) : FORTRAN 90 FORTRAN 77.,., FORTRAN 77. allocate( dd(-10:10,3) ) allocate( ee(mm) ). deallocate(dd) ; deallocate(ee),,..,..
C vs FORTRAN F:,. short int x; long int x; int x; char x; char *x; char x[7]; float x; double x; INTEGER*2 X INTEGER*4 X INTEGER*4 X BYTE X or LOGICAL*1 X CHARACTER*n X CHARACTER*7 X REAL X DOUBLE PRECISION X (), C FORTRAN.. a(3,3), a(1,1), a(2,1), a(3,1), a(1,2), a(2,2), a(3,2), a(1,3),a(2,3),a(3,3).,,, b(9)., aaa(100)1,2,3,,,100. C 0,1,2,3.99. FORTRAN bbb(0:99)., FORTRAN, bbb(1:100) bbb(100)., ccc(-100:100).
Slices, section A(:)! the whole array A(3:9)! A(m) to A(n) in steps of 1 A(3:9:1)! as above A(m:n)! A(m) to A(n) A(m:n:k)! A(m) to A(n) in steps of k A(8:3:-1)! A(8) to A(3) in steps of -1 A(8:3)! A(8) to A(3) step 1 => Zero size A(m:)! from A(m) to default UPB A(:n)! from default LWB to A(n) A(::2)! from default LWB to UPB step 2 A(m:m)! 1 element section A(m)! scalar element - not a section aaa(:,3)=bb(:) abc=maxval(aaa(:,2)) abc=minval(abs(aaa(:,:)). (?) FORTRAN 90.
program area implicit none real*8 a,b,c real*8 s,p,aaa read(5,*) a,b,c p=a+b+c s=p/2.0d0 aaa=sqrt(s*(s-a)*(s-b)*(s-c)) write(6,*) aaa stop end program area (5 ), read(5,*) (6 ), write(6,*) a.out < input > output.. implicit none real*8 a,b,c real*8 d read(5,*) a,b,c d=b**2-4.d0*a*c if(d < 0.0d0)then write(6,*) 'NO REAL ROOTS' end if if(d == 0.0d0)then write(6,*) 'two identical roots' write(6,*) -b/(2.d0*a) endif if(d> 0.0d0)then write(6,*) 'two distinct roots' write(6,*) (-b+sqrt(d))/(2.0d0*a) write(6,*) (-b-sqrt(d))/(2.0d0*a) endif stop end
ihlee@cetus0:~/lecture_programming$ a.out Enter nn 10 nn 10 55.00000000000000 aaa ihlee@cetus0:~/lecture_programming$ a.out < input > output., interactive.,.
ihlee@cetus0:~/lecture_programming$ a.out Enter nn 10 nn 10 55.00000000000000 aaa ihlee@cetus0:~/lecture_programming$,... :,,. :,. Subroutine.,.,. Subroutine.
!234567890 program i_sum implicit none integer mm real*8 bbb write(6,*) ' Enter nn' read(5,*) mm write(6,*) 'nn ',mm call int_sum(mm,bbb) write(6,*) bbb,' aaa' end program i_sum subroutine int_sum(nn,aaa) implicit none integer nn real*8 aaa integer i! write(6,*) ' Enter nn'! read(5,*) nn! write(6,*) 'nn ',nn aaa=0.0d0 do i=1,nn aaa=aaa+float(i) end do! write(6,*) aaa,' aaa' end subroutine int_sum ihlee@cetus0:~/lecture_programming$ pgf90 a.f90 sub1.f90 ihlee@cetus0:~/lecture_programming$ a.out Enter nn 100 nn 100 5050.000000000000 aaa ihlee@cetus0:~/lecture_programming$
!234567890 implicit none integer mm real*8 bbb! write(6,*) ' Enter nn' read(5,*) mm write(6,*) 'nn ',mm call int_sum(mm,bbb) write(6,*) bbb,' aaa' call inv_sum(mm,bbb) write(6,*) bbb,' aaa' end subroutine int_sum(nn,aaa) implicit none integer nn real*8 aaa integer i! write(6,*) ' Enter nn'! read(5,*) nn! write(6,*) 'nn ',nn aaa=0.0d0 do i=1,nn aaa=aaa+float(i) end do! write(6,*) aaa,' aaa' end subroutine int_sum ihlee@cetus0:~/lecture_programming$ pgf90 a.f90 sub1.f90 sub2.f90 ihlee@cetus0:~/lecture_programming$ a.out Enter nn 10 nn 10 55.00000000000000 aaa 2.928968253968254 aaa (subroutine, function).
implicit none integer n,nn real*8 a,b real(8) x,f,area,h,s integer k f(x)=x*x! write(6,*) 'a,b,n' read(5,*) a,b,n h=(b-a)/float(n) nn=n-1 s=f(a)+f(b) do k=1,nn s=s+2.0d0*f(a+float(k)*h) enddo area=h*s/2.0d0 write(6,*) a,b,area stop end ihlee@cetus0:~/lecture_programming$ a.out a,b,n 0.d0 1.d0 55 0.0000000000000000 1.000000000000000 0.3333884297520662 FORTRAN STOP ihlee@cetus0:~/lecture_programming$ ihlee@cetus0:~/lecture_programming$ a.out a,b,n 0.d0 1.d0 100 0.0000000000000000 1.000000000000000 0.3333500000000000 FORTRAN STOP ihlee@cetus0:~/lecture_programming$
pgf90 *.f90 pgf90 a.f90 b.f90 c.f90 d.f e.f pgf90 a.f90 b.o c.o d.o e.o makefile
IMPLICIT NONE? IMPLICIT NONE :.,.. REAL*8, INTEGER, CHARACTER*3, LOGICAL, COMPLEX*16. TYPO..,., USE,. IMPLICIT REAL*8 (A-H,O-Z), A-H O-Z., implicit none.... INTEGER, DIMENSION(:), ALLOCATABLE, SAVE, PRIVATE :: ivari_name.,.. INTEGER, ALLOCATABLE :: ivari_name(:) private save,. public private. implicit none private save integer... real*8... real*8, allocatable :: integer, allocatable :: logical... public ::...( ),, ().. Implicit none -.
!234567890 PROGRAM xxx_1 IMPLICIT NONE INTEGER, PARAMETER :: idp=kind(0.d0)! CHARACTER*8 fnnd ; CHARACTER*10 fnnt CHARACTER(8) fnnd ; CHARACTER(10) fnnt LOGICAL ll1,ll2 INTEGER ii,jj,i,nn INTEGER :: kk,npt! REAL*8 aa,pi REAL(8) bb,dx,xx REAL(idp) cc REAL(KIND=idp) dd REAL(KIND=idp) :: ee REAL(8) ff! COMPLEX(idp) zz1,zz2,zz3 COMPLEX(idp) :: zz4! REAL(8), ALLOCATABLE :: xvector(:),xmatrix(:,:) REAL(8), ALLOCATABLE :: yvector(:) REAL(8), ALLOCATABLE :: ymatrix(:,:) REAL(8), ALLOCATABLE :: qvector(:)! REAL(8) :: sma(3) REAL(8) smb(3)! CALL DATE_AND_TIME(DATE=fnnd,TIME=fnnt) WRITE(6,*) ' ',' date ',fnnd,' time ',fnnt,' '! INQUIRE(FILE='input_file',EXIST=ll1)! pi=4.0d0*atan(1.0d0) WRITE(6,*) 'pi ',pi pi=acos(-1.0d0) WRITE(6,*) 'pi ',pi IF(ll1)THEN WRITE(6,*) 'file is present' ELSE WRITE(6,*) 'file is absent' ENDIF WRITE(6,*) idp,' idp' ll1=.true. ll2=.false. WRITE(6,*) ll1,' ll1' WRITE(6,*) ll2,' ll2' ii=0 jj=1 kk=-1 ll1= (1 == 2) ll2= (1 /= 2) IF(ii >= jj) ll1=.true. IF(ii == jj) ll2=.false. WRITE(6,*) ii,' ii' WRITE(6,*) jj,' jj' aa=0.0_idp! note bb=0.0d0 cc=0.0d0 dd=0.0d0 zz1=0.0_idp zz2=cmplx(1.0_idp,2.0_idp)! note zz3=zz1+zz2 zz3=zz1-zz2 zz3=zz1*zz2 zz3=zz1/zz2 zz3=zz1**2 aa=real(zz3) bb=imag(zz3) cc=abs(zz3) zz4=conjg(zz3) WRITE(6,*) aa,' aa' WRITE(6,*) bb,' bb' WRITE(6,*) cc,' cc' WRITE(6,*) dd,' dd' WRITE(6,*) zz1,' zz1' WRITE(6,*) zz2,' zz2' OPEN(1,FILE='output1',FORM='FORMATTED') WRITE(1,'(i8)') ii WRITE(1,'(2i8)') ii,jj WRITE(1,'(i8,2x,i8)') ii,jj WRITE(1,'(i8,2x,i8,3x,a5)') ii,jj,'ii,jj' WRITE(1,*) aa,' aa' WRITE(1,*) bb,' bb' WRITE(1,*) cc,' cc' WRITE(1,*) dd,' dd' WRITE(1,*) zz1,' zz1' WRITE(1,*) zz2,' zz2' CLOSE(1)
OPEN(2,FILE='output2',FORM='FORMATTED') WRITE(2,'(i8)') ii WRITE(2,'(2i8)') ii,jj CLOSE(2) sma=0.0d0 smb=1.0d0 ii=10 jj=10 nn=ii ALLOCATE(xvector(nn)) ; ALLOCATE(yvector(nn)) ALLOCATE(xmatrix(nn,nn)) ALLOCATE(ymatrix(nn,nn)) ALLOCATE(qvector(0:nn)) IF(ALLOCATED(xvector))THEN WRITE(6,*) 'xvector is allocated' WRITE(6,*) 'size ', SIZE(xvector) DO i=1,nn xvector(i)=float(i) ENDDO yvector=-xvector ENDIF WRITE(6,*) 'sum ', SUM(xvector) WRITE(6,*) 'sum ', SUM(yvector) yvector=sqrt(abs(xvector)) IF(ALLOCATED(xmatrix))THEN WRITE(6,*) 'xmatrix is allocated' WRITE(6,*) 'size ', SIZE(xmatrix) ENDIF xvector=0.0d0 yvector=0.0d0 xmatrix=0.0d0 ymatrix=transpose(xmatrix) ymatrix(:,1)=yvector(:) yvector=matmul(xmatrix,xvector) dd=dot_product(xvector,yvector) dd=maxval(xvector) ee=minval(yvector) dd=maxval(abs(xvector)) dd=maxval(abs(xvector-yvector))/(maxval(abs(xvector))+maxval(abs(yvector))+1.0d-8) aa=sum(xvector)/float(size(xvector)) bb=product(yvector) DEALLOCATE(xvector) ; DEALLOCATE(yvector) DEALLOCATE(xmatrix) DEALLOCATE(ymatrix) DEALLOCATE(qvector) INQUIRE(FILE='del',EXIST=ll1) IF(LL1)THEN WRITE(6,*) 'del is present' OPEN(11,FILE='del') CLOSE(31,STATUS='DELETE') ENDIF aa=0.0d0 bb=1.0d0 npt=100000000 dx=(bb-aa)/float(npt-1) cc=0.0d0 DO i=1,npt xx=aa+float(i-1)*dx cc=cc+ff(xx) ENDDO cc=cc*dx write(6,*) cc,' cc' STOP END PROGRAM xxx_1 FUNCTION ff(x) REAL(8) ff REAL(8) x ff=x*x RETURN END FUNCTION ff
ihlee@cetus0:~/lecture_programming$ a.out date 20080109 time 152244.148 pi 3.141592653589793 pi 3.141592653589793 file is absent 8 idp T ll1 F ll2 0 ii 1 jj 0.0000000000000000 aa 0.0000000000000000 bb 0.0000000000000000 cc 0.0000000000000000 dd (0.0000000000000000,0.0000000000000000) zz1 (1.000000000000000,2.000000000000000) zz2 xvector is allocated size 10 sum 55.00000000000000 sum -55.00000000000000 xmatrix is allocated size 100 0.3333333283332683 cc FORTRAN STOP ihlee@cetus0:~/lecture_programming$
20 TOP 10 http://incredible.egloos.com/3216474.,.,.,. Fortran Calls C http://docs.sun.com/app/docs/doc/801-5492/6hvudcqa3?l=ko&a=view C Calls Fortran http://docs.sun.com/app/docs/doc/801-5492/6hvudcqa4?l=ko&a=view
Fortran Calls C void simref_ ( b4, i4, r4, d8 ) /* Simple types, passed by reference, from f90 (f90 calls C)*/ int * b4 ; int * i4 ; float * r4 ; double * d8 ; { *b4 = 1 ; *i4 = 9 ; *r4 = 9.9f ; *d8 = 9.9F ; } PROGRAM SimpleRef! Pass some simple types, by reference, to C (f90 calls C) LOGICAL b4! Default kind is 4-byte INTEGER i4! Default kind is 4-byte REAL r4! Default kind is 4-byte DOUBLE PRECISION d8! This is 8-byte CALL SimRef ( b4, i4, r4, d8 ) WRITE( *, '(L4,I4,F6.1,F6.1)') b4, i4, r4, d8 END PROGRAM SimpleRef ihlee@cetus0:~/lecture_programming$ vi simref.c ihlee@cetus0:~/lecture_programming$ vi simplref.f90 ihlee@cetus0:~/lecture_programming$ cc -c simref.c ihlee@cetus0:~/lecture_programming$ pgf90 simplref.f90 simref.o ihlee@cetus0:~/lecture_programming$ a.out T 9 9.9 9.9 http://docs.sun.com/app/docs/doc/801-5492/6hvudcqa3?l=ko&a=view
,,.,..,.??. : (class) : (object, instance)
-C ( ) -c object. real*8 a(3,2),yone(6) 3*2=6. (row) (column). 1. a(1,1) a(2,1) a(3,1) a(1,2) a(2,2) a(3,2)ij=3*(j-1)+i a(i,j) (yone,, 6.). do i=1,3 do j=1,2 ij=3*(j-1)+i yone(ij)=a(i,j) enddo enddo akvector(natom,3,nk) 1 ( ). do iatom=1,natom kp=3*nk*(iatom-1) do k=1,nk xone(kp+3*(k-1)+1)=akvector(iatom,1,k) xone(kp+3*(k-1)+2)=akvector(iatom,2,k) xone(kp+3*(k-1)+3)=akvector(iatom,3,k) enddo enddo..,, 1..
Subroutine, function( ),. double precision, double precision.. (.) A(3,2). subroutine abc1(a) implicit none real*8 a(3,1) subroutine abc2(a) implicit none real*8 a(3,2).,.,.,. (.) subroutine abc3(a) implicit none real*8 a(1,2) ( 77) subroutine abc4(a,lda) implicit none integer lda real*8 a(lda,1). a(3,3), a(1,1), a(2,1), a(3,1), a(1,2), a(2,2), a(3,2), a(1,3),a(2,3),a(3,3).,,, b(9).
implicit none integer m integer j! read(5,*) m do j=2,m call test_prime(j) end do stop end subroutine test_prime(n) implicit none integer n integer k if( n <=1 ) return if( n == 2) write(6,*) n,' is a prime' if( n == 2) return k=2 30 if (n/k*k == n) goto 70 k=k+1 if(k <= n/2) go to 30 write(6,*) n,' is a prime' return 70 continue write(6,*) n, ' is not a prime' return end ihlee@cetus0:~/lecture_programming$ a.out 20 2 is a prime 3 is a prime 4 is not a prime 5 is a prime 6 is not a prime 7 is a prime 8 is not a prime 9 is not a prime 10 is not a prime 11 is a prime 12 is not a prime 13 is a prime 14 is not a prime 15 is not a prime 16 is not a prime 17 is a prime 18 is not a prime 19 is a prime 20 is not a prime FORTRAN STOP ihlee@cetus0:~/lecture_programming$
!234567890 implicit none integer nm,n,matz,ierr real*8, allocatable :: ar(:,:),ai(:,:),w(:),zr(:,:),zi(:,:),fv1(:),fv2(:),fm1(:,:) real*8 tmp integer i,j,k matz=0 matz=1 nm=3 n=nm open(1,file='input',form='formatted') read(1,*) n close(1) write(6,*) n,' n' allocate(ar(nm,n),ai(nm,n),w(n),zr(nm,n),zi(nm,n),fv1(n),fv2(n),fm1(2,n)) do i=1,n do j=1,i! lower-triangular-form.!------[ set up the hermitian matrix : lower triangular form ar(i,j)=dfloat(i+j)**2 ai(i,j)=1.d0/dfloat(i-j) if(i == j) ai(i,j)=0.0d0 hermitian enddo enddo do i=1,n do j=1,i ar(j,i)=ar(i,j) ai(j,i)=-ai(i,j) enddo enddo call ch(nm,n,ar,ai,w,matz,zr,zi,fv1,fv2,fm1,ierr) write(6,*) 'eigenvalues' do i=1,n write(6,*) w(i) enddo write(6,*) 'eigenvectors' do j=1,n write(6,*) (dcmplx(zr(k,j),zi(k,j)),k=1,n) enddo do i=1,n do j=1,n tmp=0.0d0 do k=1,n tmp=tmp+dcmplx(zr(k,i),-zi(k,i))*dcmplx(zr(k,j),zi(k,j)) enddo write(6,*) i,j,tmp enddo endd deallocate(ar,ai,w,zr,zi,fv1,fv2,fm1)!. stop end 3 3 n eigenvalues -3.295706627065439 0.2569711761519913 59.03873545091345 eigenvectors (0.7426815277007046,0.2134499083118499) (0.3321344158486125,-0.1647469485075425) (-0.5151780835984430,0.000000000000000) (-0.4017627191549581,-0.3764084131040000) (0.7336179945919501,0.2195919189859811) (-0.3323965270105607,0.000000000000000) (0.3152758193792096,-1.9179854283106629E-002) (0.5252646650262870,-1.5040996247447431E-002) (0.7900025892433221,0.000000000000000) 1 1 0.9999999999999998 1 2-2.6099456797157305E-016 1 3-2.2218013019659200E-016 2 1-2.6099456797157305E-016 2 2 0.9999999999999998 2 3-2.0545631168600309E-016 3 1-2.2218013019659200E-016 3 2-2.0545631168600309E-016 3 3 0.9999999999999996
program conjugate_gradient implicit none integer iters, its, n logical :: converged real*8 :: tol, up, alpha, beta real*8, allocatable :: a(:,:), b(:), x(:), r(:), u(:), p(:), xnew(:) read(5,*) n, tol, its allocate( a(n,n), b(n), x(n), r(n), u(n), p(n), xnew(n) ) open(10, file='data') read(10,*) a read(10,*) b close(10) x=1.d0! 1.d0. r=b-matmul(a,x)!.. p=r!. iters=0 do! do loop. iters=iters+1 u=matmul(a,p) up=dot_product(r,r)! (inner product). alpha=up/dot_product(p,u) xnew=x+p*alpha r=r-u*alpha!. beta=dot_product(r,r)/up p=r+p*beta! beta. converged=(maxval(abs(xnew-x))/maxval(abs(x)) < tol)!. x=xnew!. if( converged.or. iters == its) exit! do. enddo write(6,*) iters write(6,*) x deallocate( a,b,x,r,u,p,xnew ) end program conjugate_gradient!.
module common_data implicit none private!.. ( default.) save integer iii1,jjj1 real*, allocatable :: aa1(:,:),bb1(:) logical first_call!.!..true.,.false.. character*2, allocatable :: ccc1(:)... public :: aa1,bb1,iii1!,. (.) end module common_data,.,,,, public. public.,, default public. module common_data2 implicit none private save... public ::... end module common_data2 :. :. Interface.,,. : http://incredible.egloos.com/2950367
( ) CONTAINS (, ). END FUNCTION SUBROUTINE. End function fname,, end subroutine sname... http://www.liv.ac.uk/hpc/htmlhpfcourse/htmlhpfcourseslidesnode142.html
..!234567890 MODULE wf_wr IMPLICIT NONE save public integer nsim real*8 temper real*8, allocatable :: wf(:),wr(:) END MODULE wf_wr!. USE wf_wr IMPLICIT NONE real*8 tol,sol,sol1,sol2 real*8 tmp integer i real*8 bennetteq,zeroin temper=1.d0 ; nsim=10 allocate(wf(nsim),wr(nsim)) wf=1. ; wr=-1. do i=1,nsim call random_number(tmp) wf(i)=wf(i)+0.2*(tmp-0.5) call random_number(tmp) wr(i)=wr(i)+0.2*(tmp-0.5) enddo sol1=maxval(wf) ; sol2=minval(wf) write(6,*) sol1,bennetteq(sol1) write(6,*) sol2,bennetteq(sol2) tol=1.d-16 sol=zeroin(sol1,sol2,bennetteq,tol)! netlib. root finding. : bennetteq. write(6,*) sol,bennetteq(sol) deallocate(wf,wr) stop end
USE, ONLY :, public,. USE USE ONLY.. save, public ( private),. implicit none... 90. implicit none.,,,,,.. -C.. USE stack, local_use_pop => pop pop, stack pop local_pop. USE stack, ONLY: pos, local_pop=> pop stack pos pop., pop local_pop.. makefile., ximage_lbfgs.f90 timpose.f90.. USE aaa USE bbb! bbb public,, (function, subroutine). USE ccc USE ddd, ONLY: qqq,ppp! ddd,, (subroutine, function).
implicit none,, ()..,...,., = {}+{ }. (contain)..,. (, ) USE. 77,, common common. 90 77...,.,.,.,,,, public. public.,, default public.
cmpl=pgf90 segl=pgf90 Ocft= -c -fast cft= -c #LINK=-L/usr/lib/xmm -L/usr/lib/xmm/atlas LINK= -L/usr/lib/sse2 LIBS = -llapack -lf77blas -latlas -lf2c AMD= amd.o predictor_corrector_nose.o energy_force.o tight_binding_lapack_c.o \ mini_driver.o minimization.o lbfgs.o bond_boost.o erf.o maxwell_boltzmann.o \ check_transition.o rmarin.o angle_boost.o dihedral.o dihedral_boost.o amd.x: $(AMD) $(segl) -o amd.x $(AMD) $(LINK) $(LIBS) amd.o:amd.f90 predictor_corrector_nose.o bond_boost.o angle_boost.o dihedral_boost.o $(cmpl) $(Ocft) amd.f90 energy_force.o:energy_force.f90 tight_binding_lapack_c.o $(cmpl) $(Ocft) energy_force.f90 tight_binding_lapack_c.o:tight_binding_lapack_c.f90 $(cmpl) $(Ocft) tight_binding_lapack_c.f90 predictor_corrector_nose.o:predictor_corrector_nose.f90 $(cmpl) $(Ocft) predictor_corrector_nose.f90 mini_driver.o:mini_driver.f90 $(cmpl) $(Ocft) mini_driver.f90 minimization.o:minimization.f90 $(cmpl) $(Ocft) minimization.f90 bond_boost.o:bond_boost.f90 $(cmpl) $(Ocft) bond_boost.f90 angle_boost.o:angle_boost.f90 $(cmpl) $(Ocft) angle_boost.f90 dihedral_boost.o:dihedral_boost.f90 dihedral.o $(cmpl) $(Ocft) dihedral_boost.f90 dihedral.o:dihedral.f90 erf.o:erf.f $(cmpl) $(Ocft) erf.f clean: rm -f *.x *.o *.mod *.M core* touch: touch *.f90 *.i makefile ; chmod 600 *.f90 *.i makefile ; ls -l *.f90 *.i makefile rmo: rm -f *.o *.mod *.M core* lsl: ls -l *.f90 makefile *.i, makefile. make. ls -ltra. ls -Fs make clean make
! module metric implicit none real*8 aa,bb,cc,dd end module metric program abc USE metric implicit none real*8 sol1,sol2,tol real*8 sol,x real*8, external :: func,zeroin integer i aa=0.d0 bb=-2.d0 cc= 1.d0 dd=2.d0 sol1=-1.0d0 sol2= 0.1d0 write(6,*) sol1,func(sol1) write(6,*) sol2,func(sol2) tol=1.d-16 ; sol=zeroin(sol1,sol2,func,tol) write(6,*) sol,func(sol),' sol, f(sol)' sol1=-3.0d0 sol2=-1.9d0 write(6,*) sol1,func(sol1) write(6,*) sol2,func(sol2) tol=1.d-16 ; sol=zeroin(sol1,sol2,func,tol) write(6,*) sol,func(sol),' sol, f(sol)' sol1= 0.9d0 sol2= 1.1d0 write(6,*) sol1,func(sol1) write(6,*) sol2,func(sol2) tol=1.d-16 ; sol=zeroin(sol1,sol2,func,tol) write(6,*) sol,func(sol),' sol, f(sol)' sol1= 1.9d0 sol2= 2.1d0 write(6,*) sol1,func(sol1) write(6,*) sol2,func(sol2) tol=1.d-16 ; sol=zeroin(sol1,sol2,func,tol) write(6,*) sol,func(sol),' sol, f(sol)' stop end program abc real*8 function func(x) USE metric implicit none real*8 x func=(x-dd)*(x-aa)*(x-bb)*(x-cc) return end -1.000000000000000-6.000000000000000 0.1000000000000000 0.3591000000000000 2.7076140240079107E-019 1.0830456096031642E-018 sol, f(sol) -3.000000000000000 60.00000000000000-1.900000000000000-2.148900000000002-2.000000000000000 0.0000000000000000 sol, f(sol) 0.9000000000000000 0.2870999999999999 1.100000000000000-0.3069000000000003 1.000000000000000 0.0000000000000000 sol, f(sol) 1.900000000000000-0.6669000000000005 2.100000000000000 0.9471000000000009 2.000000000000000 0.0000000000000000 sol, f(sol) FORTRAN STOP.,,.,, (.). Cf. 77 common. : http://incredible.egloos.com/3040113
module dihedral_angle! Written by In-Ho Lee, KRISS, October 4 (2003). implicit none private save real*8 diangle,dadi(3),dadj(3),dadk(3),dadl(3) public :: mk_diangle,diangle,dadi,dadj,dadk,dadl contains subroutine mk_diangle(ri,rj,rk,rl) implicit none real*8 ri(3),rj(3),rk(3),rl(3) real*8 rij(3),rkj(3),rkl(3),xnorm_ij_kj,xnorm_kj_kl real*8 rij_kj(3),rkj_kl(3),rtmp(3),xnorm_kj,xnorm_mj,xnorm_nk real*8 rmj(3),rnk(3),xxr,xxi real*8 argu,qsign rij=ri-rj ; rkj=rk-rj ; rkl=rk-rl call cross(rij,rkj,rij_kj) call cross(rkj,rkl,rkj_kl) call cross(rij_kj,rkj_kl,rtmp) xnorm_ij_kj=sqrt(dot_product(rij_kj,rij_kj)) xnorm_kj_kl=sqrt(dot_product(rkj_kl,rkj_kl)) argu=dot_product(rij_kj,rkj_kl)/(xnorm_ij_kj*xnorm_kj_kl) ; argu=min(1.0d0,max(-1.0d0,argu)) qsign=1.0d0 ; if(dot_product(rkj,rtmp) < 0.0d0) qsign=-1.0d0 diangle=qsign*acos(argu) call cross(rij,rkj,rmj) call cross(rkj,rkl,rnk) xnorm_kj=sqrt(dot_product(rkj,rkj)) xnorm_mj=sqrt(dot_product(rmj,rmj)) xnorm_nk=sqrt(dot_product(rnk,rnk)) dadi(:)=rmj(:)*xnorm_kj/xnorm_mj**2 dadl(:)=-rnk(:)*xnorm_kj/xnorm_nk**2 xxr=dot_product(rij,rkj)/xnorm_kj**2-1.0d0 xxi=dot_product(rkl,rkj)/xnorm_kj**2 dadj(:)=xxr*dadi(:)-xxi*dadl(:) xxr=dot_product(rkl,rkj)/xnorm_kj**2-1.0d0 xxi=dot_product(rij,rkj)/xnorm_kj**2 dadk(:)=xxr*dadl(:)-xxi*dadi(:) dadi=dadi*qsign dadj=dadj*qsign dadk=dadk*qsign dadl=dadl*qsign end subroutine mk_diangle subroutine cross(r1,r2,r3) implicit none real*8 r1(3),r2(3),r3(3) r3(1) = r1(2)*r2(3)-r1(3)*r2(2) r3(2) = r1(3)*r2(1)-r1(1)*r2(3) r3(3) = r1(1)*r2(2)-r1(2)*r2(1) end subroutine cross end module dihedral_angle Dihedral angle
!234567890 implicit none real*8, external :: func real*8 rslt,aa,bb integer n n=100 aa=0.0d0 bb=1.0d0 call simpson(func,n,aa,bb,rslt) write(6,*) rslt,' rslt' stop end subroutine simpson(func,n,aa,bb,rslt) implicit none real*8 func integer n real*8 rslt,aa,bb real*8 h,xx integer j logical lodd h=(bb-aa)/float(n) rslt=(func(aa)+func(bb)) lodd=.true. do j=1,n-1 xx=aa+h*float(j) if(lodd)then rslt=rslt+4.0d0*func(xx) else rslt=rslt+2.0d0*func(xx) endif lodd=(.not. lodd) enddo rslt=rslt*h/3.0d0 return end real*8 function func(x) implicit none real*8 x func=x*x return end ---------------------------------------------- 0.3333333333333334 rslt FORTRAN STOP n even number. Composite Simpson's rule : [a,b] step length: h=(b-a)/n x_i=a+i h i=0,1,2,3,...,n-1,n n. x_0=a x_n=b : http://incredible.egloos.com/3021478
, Fornberg ( ) Uniform grid. (nonuniform grid),. (.),. 0 =(k=0) k=1 () k grid. f(x) : http://incredible.egloos.com/3032905
subroutine weights(xi,x,n,m,c)! purpose: computes weights for the mth derivative on arbitrarily! spaced points! Jan 04, 1996 declare all variables! date: Nov 15, 1995 (copied from Fornberg)! written by In-Ho Lee, Beckman Inst., UIUC, April 13 1997.! +real_constant=double option in HP f90! +--------------------------------------------------------------------+! INPUT PARAMETERS:! xi point at which the approximations are to be accurate! x x-coordinates for the grid points, array dimensioned x(0:n)! n the grid points are at x(0),x(1),...x(n) (i.e. n+1 in all)! m highest order of derivative to be approximated!! OUTPUT PARAMETERS! c weights, array dimensioned c(0:n,0:n,0:m)! on return, the element c(j,k,i) contains the weight to be! applied at x(k) when the i:th derivative is approximated! by a stencil extending over x(0),x(1),...,x(j).! +--------------------------------------------------------------------+! implicit real*8 (a-h,o-z) dimension x(0:n),c(0:n,0:n,0:m) c(0,0,0)= 1.0d0 c1 = 1.0d0 c4 = x(0)-xi do 40 j=1,n mn = min(j,m) c2 = 1.0d0 c5 = c4 c4 = x(j)-xi do 20 k=0,j-1 c3=x(j)-x(k) c2=c2*c3 if (j.le.m) c(k,j-1,j)=0.0d0 c(k,j,0) =c4*c(k,j-1,0)/c3 do 10 i=1,mn 10 c(k,j,i)=(c4*c(k,j-1,i)-i*c(k,j-1,i-1))/c3 20 continue c(j,j,0)=-c1*c5*c(j-1,j-1,0)/c2 do 30 i=1,mn 30 c(j,j,i)=c1*(i*c(j-1,j-1,i-1)-c5*c(j-1,j-1,i))/c2 40 c1=c2 return end, Fornberg ( ) Fornberg Finite difference., high-order.. Finite-difference. 10.,. (, : ).,, m.,,... Finite difference: http://en.wikipedia.org/wiki/finite_difference http://mathworld.wolfram.com/finitedifference.html http://en.wikipedia.org/wiki/numerical_differentiation
!234567890 module rk_df implicit none public save integer neq,info(15),idid,lrw,liw real*8 t, tout integer, allocatable :: iwork(:),ipar(:) real*8, allocatable :: y(:), rwork(:), rtol(:), atol(:), rpar(:) contains subroutine initial implicit none t=0.0d0 tout=0.1d0 neq=3 allocate(y(neq)) y(1)=0.d0 y(2)=1.d0 y(3)=1.d0 info(1)=0 lrw= 33+7*neq ; liw= 34 allocate(iwork(liw)) ; allocate(rwork(lrw)) allocate(atol(neq),rtol(neq)) rtol=1.d-10 ; atol=1.d-10 allocate(rpar(neq)) ; allocate(ipar(neq)) end subroutine initial subroutine final implicit none deallocate(y) ; deallocate(rwork) ; deallocate(iwork) deallocate(atol,rtol) ; deallocate(rpar) ; deallocate(ipar) end subroutine final! DF(X,U,UPRIME,RPAR,IPAR) x, u.! uprime(). subroutine df(x,u,uprime,rpar_z,ipar_z) implicit none real*8 x,u(1) real*8 rpar_z(1) integer ipar_z(1) real*8 uprime(1) uprime(1)=u(2)*u(3) uprime(2)=-u(1)*u(3) uprime(3)=-0.522d0*u(1)*u(2) end subroutine df end module rk_df! program rk_test USE rk_df, ONLY : neq, t, y, tout, info, rtol, atol, idid, rwork, lrw, iwork, liw, rpar, ipar USE rk_df, ONLY : df,initial,final implicit none call initial call dderkf(df, neq, t, y, tout, info, rtol, atol, idid, rwork, lrw, iwork, liw, rpar, ipar)! write(6,*) t,' t' write(6,*) tout,' tout' write(6,*) y,' y' write(6,'(i5,1x,a4)') idid,'idid' write(6,'(15i5,1x,a4)') info,'info' tout=tout+0.1d0 call dderkf(df, neq, t, y, tout, info, rtol, atol, idid, rwork, lrw, iwork, liw, rpar, ipar)! write(6,*) t,' t' write(6,*) tout,' tout' write(6,*) y,' y' write(6,'(i5,1x,a4)') idid,'idid' write(6,'(15i5,1x,a4)') info,'info' tout=tout+0.1d0 call dderkf(df, neq, t, y, tout, info, rtol, atol, idid, rwork, lrw, iwork, liw, rpar, ipar)! write(6,*) t,' t' write(6,*) tout,' tout' write(6,*) y,' y' write(6,'(i5,1x,a4)') idid,'idid' write(6,'(15i5,1x,a4)') info,'info' call final end program rk_test allocate(vector(-89:100)). -89, -88, -87,...98, 99, 100., allocate(vector(100)) 100 vector., 1 100. allocate(vector(1:100)). 1. : http://incredible.egloos.com/3143433
implicit none integer j,i integer narr,mref real*8, allocatable :: xarr(:,:) real*8, allocatable :: yarr(:,:) logical lexist narr=10 mref=3 allocate(xarr(narr,mref)) allocate(yarr(narr,mref)) inquire(file='fort.2',exist=lexist) if(lexist)then open(2,file='fort.2',access='direct',recl=8*narr,form='unformatted') do j=1,mref read(2,rec=j) yarr(:,j) enddo close(2) endif do j=1,mref do i=1,narr xarr(i,j)=float(i)**2/float(j) enddo enddo write(6,*) maxval(abs(xarr-yarr)) open(2,file='fort.2',access='direct',recl=8*narr,form='unformatted') do j=1,mref write(2,rec=j) xarr(:,j) enddo close(2) deallocate(xarr) deallocate(yarr) stop end implicit none integer j,i integer narr,mref real*8, allocatable :: xarr(:,:) real*8, allocatable :: yarr(:,:) logical lexist narr=10 mref=3 allocate(xarr(narr,mref)) allocate(yarr(narr,mref)) inquire(file='fort.12',exist=lexist) if(lexist)then open(12,file='fort.12',form='unformatted') read(12) yarr(:,:) close(12) endif do j=1,mref do i=1,narr xarr(i,j)=float(i)**2/float(j) enddo enddo write(6,*) maxval(abs(xarr-yarr)) open(12,file='fort.12',form='unformatted') write(12) xarr(:,:) close(12) deallocate(xarr) deallocate(yarr) stop end : http://incredible.egloos.com/3044815
READ(*,100) i,j WRITE(*,100) i,j READ(*,FMT=200) x,y WRITE(*,200) x,y 100 FORMAT(2I) 200 FORMAT(2F10.6) READ(*,'(2I)') i,j WRITE(*,'(2F12.6)') x,y WRITE(6,`(2F12.6)`) 12.6, -131.45678 INTEGER :: i, j REAL, DIMENSION(10) :: a REAL, DIMENSION(10,10) :: b READ(*,*) (a(j),j=1,10) WRITE(*,*) (a(j), j=10,1,-1) WRITE(*,*) ((b(i,j),i=1,10), j=1,10). formatted, unformatted. Unformatted. /... form= formatted serial, direct. Direct CD, 3 7. access= direct,. recl. Record length.
open(11,file= abcd.txt, form= unformatted, access= direct, recl=8*nnsszz) Unformatted... write(11,rec=k) (aa(i),i=1,nnsszz) Unit number... read(11,rec=kk) (bb(j),j=1,nnsszz) /(record length). Unformatted read/write... cf. netcdf, HDF hdf.ncsa.uiuc.edu/products/hd5/index.html
Inquiry list EXIST=lex!.true. or.false. OPENED=lod!.true. or.false. NUMBER=unum! unit number NAME=fnm! filename ACCESS=acc! 'DIRECT' or 'SEQUENTIAL' SEQUENTIAL=seq! 'YES or 'NO' DIRECT=dir! 'YES' or 'NO' FORMATTED=fmt! 'YES' or 'NO' UNFORMATTED=unfmt! 'YES' or 'NO' FORM=frm! 'FORMATTED' or 'UNFORMATTED' NEXTREC=recn! number of next record RECL=recl! record length inquire(file='fort.40',exist=lexist) if(lexist)then else. endif allocated(abc) abc..true.,.false.. if(allocated(abc)) deallocate(abc) if(.not. allocated(abc)) allocate(abc(0:100)).
implicit none real*8, allocatable :: xdata(:),ydata(:) character*15, file_a,file_b integer j,jtemp,ndata real*8 rr,prob,zz write(6,*) 'Enter: ndata' read(5,*) ndata write(6,*) 'Enter: file_a' read(5,*) file_a write(6,*) file_a write(6,*) 'Enter: file_b' read(5,*) file_b write(6,*) file_b allocate(xdata(ndata),ydata(ndata)) open(11,file=file_a,form='formatted') do j=1,ndata read(11,*) jtemp,xdata(j) enddo close(11) open(21,file=file_b,form='formatted') do j=1,ndata read(21,*) jtemp,ydata(j) enddo close(21) call pearsn(xdata,ydata,ndata,rr,prob,zz) write(6,*) rr,rr**2, ' r, r^2' write(6,*) prob, ' prob' write(6,*) zz, ' z' deallocate(xdata,ydata) STOP END Linear correlation -1.1 : http://incredible.egloos.com/1992323
implicit none integer j, k,i real*8 tmp real cpu_diff real cpu_time1 real cpu_time2 real ( kind = 8 ) real_diff real ( kind = 8 ) real_time1 real ( kind = 8 ) real_time2 call cpu_time (cpu_time1 )! call real_time ( real_time1 ) tmp=0.0d0 do i=1,1000 do j=1,1000 do k=1,100 tmp=tmp+sin(float(j*i))*cos(float(i*j)) enddo enddo enddo ihlee@cetus0:~/lecture_programming$ a.out 16.00000 FORTRAN STOP ihlee@cetus0:~/lecture_programming$ ihlee@cetus0:~/lecture_programming$ time a.out 16.00000 FORTRAN STOP real 0m16.031s user 0m16.028s sys 0m0.001s ihlee@cetus0:~/lecture_programming$ call cpu_time ( cpu_time2 )! call real_time ( real_time2 ) cpu_diff = cpu_time2 - cpu_time1! real_diff = real_time2 - real_time1 write(6,*) cpu_diff! write(6,*) real_diff stop end : http://incredible.egloos.com/3047318
subroutine timestamp ( )!*****************************************************************************80!!! TIMESTAMP prints the current YMDHMS date as a time stamp.!!example:!! May 31 2001 9:45:54.872 AM!! Modified:!! 31 May 2001!! Author:!! John Burkardt!! Parameters:!! None! implicit none character ( len = 8 ) ampm integer d character ( len = 8 ) date integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s character ( len = 10 ) time integer values(8) integer y character ( len = 5 ) zone call date_and_time ( date, time, zone, values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) 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!234567890 implicit none call timestamp() stop end ----------------------------------------- March 12 2007 1:58:40.868 PM FORTRAN STOP write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end : http://incredible.egloos.com/3047077
subroutine lbfgs_mini(object,posi,grad,nn)! Written by In-Ho Lee, KRISS, March 3 (2003). implicit none integer nn real*8 object,posi(nn),grad(nn) integer n_local,m,msave real*8, allocatable :: diag(:),w(:) real*8 eps,xtol,gtol,t1,t2,stpmin,stpmax integer iprint(2),iflag,icall,mp,lp,nwork integer jnatom,j logical diagco real*8, allocatable :: xyz(:,:),fxyz(:,:)! The driver for LBFGS must always declare LB2 as EXTERNAL external lb2 common /lb3/mp,lp,gtol,stpmin,stpmax n_local=nn m=5 ; msave=7 ; nwork=n_local*(2*msave+1)+2*msave allocate(diag(n_local),w(nwork)) iprint(1)= 1 ; iprint(2)= 0!! We do not wish to provide the diagonal matrices Hk0, and! therefore set DIAGCO to FALSE. diagco=.false. ; eps= 1.0d-4 ; xtol= 1.0d-16 ; icall=0 ; iflag=0! jnatom=nn/3 allocate(xyz(jnatom,3),fxyz(jnatom,3))! 20 continue do j=1,jnatom xyz(j,1)=posi(3*(j-1)+1) xyz(j,2)=posi(3*(j-1)+2) xyz(j,3)=posi(3*(j-1)+3) enddo call energy_force(xyz,fxyz,object,jnatom) do j=1,jnatom grad(3*(j-1)+1)=-fxyz(j,1) grad(3*(j-1)+2)=-fxyz(j,2) grad(3*(j-1)+3)=-fxyz(j,3) enddo! call lbfgs(n_local,m,posi,object,grad,diagco,diag,iprint,eps,xtol,w,iflag) if(iflag <= 0) go to 50 icall=icall + 1! We allow at most 2000 evaluations of Function and Gradient if(icall > 2000) go to 50 go to 20 50 continue deallocate(diag,w) deallocate(xyz,fxyz) return end ( ) 2 1 lbfgs 1. 2. Gradient, function : : http://incredible.egloos.com/3042973
!234567890 implicit none integer nm,n,matz,ierr real*8, allocatable :: ar(:,:),ai(:,:),w(:),zr(:,:),zi(:,:),fv1(:),fv2(:),fm1(:,:) real*8 tmp integer i,j,k matz=0 matz=1 nm=3 n=nm open(1,file='input',form='formatted') read(1,*) n close(1) write(6,*) n,' n' allocate(ar(nm,n),ai(nm,n),w(n),zr(nm,n),zi(nm,n),fv1(n),fv2(n),fm1(2,n)) do i=1,n do j=1,i! lower-triangular-form.!------[ set up the hermitian matrix : lower triangular form ar(i,j)=dfloat(i+j)**2 ai(i,j)=1.d0/dfloat(i-j)!------] if(i == j) ai(i,j)=0.0d0 hermitian enddo enddo do i=1,n do j=1,i ar(j,i)=ar(i,j) ai(j,i)=-ai(i,j) enddo enddo call ch(nm,n,ar,ai,w,matz,zr,zi,fv1,fv2,fm1,ierr) write(6,*) 'eigenvalues' do i=1,n write(6,*) w(i) enddo write(6,*) 'eigenvectors' do j=1,n write(6,*) (dcmplx(zr(k,j),zi(k,j)),k=1,n) enddo do i=1,n do j=1,n tmp=0.0d0 do k=1,n tmp=tmp+dcmplx(zr(k,i),-zi(k,i))*dcmplx(zr(k,j),zi(k,j)) enddo write(6,*) i,j,tmp enddo enddo deallocate(ar,ai,w,zr,zi,fv1,fv2,fm1)!. stop end Matrix diagonalization( )
!234567890 implicit none integer i,j,k character*32 filename character*4 fn integer nsize nsize=4 do j=0,3 call numeral (j,fn,nsize) write(6,*) fn filename=' del'//fn write(6,*) filename! call system("cp admd.seq"//filename) end do stop end ihlee@cetus0:~/lecture_programming$ a.out 0000 del0000 0001 del0001 0002 del0002 0003 del0003 FORTRAN STOP ihlee@cetus0:~/lecture_programming$ : http://incredible.egloos.com/3474225
!234567890 implicit none integer i,j,k character*32 subcommand character*32 subcommand2 character*4 fn character*9 fm integer nsize nsize=4 do j=0,3 call numeral (j,fn,nsize) write(6,*) fn fm=' admd'//fn subcommand=fm//'.pdb' write(6,*) subcommand! call system("cp admd.pdb"//subcommand) subcommand2=fm//'.dssp' subcommand2=trim(subcommand)//fm//'.dssp' write(6,*) subcommand2! call system("~/dssp/dsspcmbi "//subcommand2) end do stop end : http://incredible.egloos.com/3474225
MODULE cartesian TYPE point REAL :: x, y END TYPE point INTERFACE OPERATOR (.DIST. ) MODULE PROCEDURE dist END INTERFACE CONTAINS REAL FUNCTION dist( a, b ) TYPE(point), INTENT(IN) :: a, b dist = SQRT( (a%x-b%x)**2 + (a%y-b%y)**2 ) END FUNCTION dist END MODULE cartesian ihlee@cetus0:~/lecture_programming$ a.out 1.414214 FORTRAN STOP ihlee@cetus0:~/lecture_programming$ program over2 USE cartesian IMPLICIT NONE TYPE(point) :: a, b REAL distance a%x=0.0d0 a%y=0.0d0 b%x=1.0d0 b%y=1.0d0 distance = a.dist. b write(6,*) distance stop end program over2
MODULE cartesian TYPE point REAL*8 :: x, y, z END TYPE point INTERFACE OPERATOR (.DIST. ) MODULE PROCEDURE dist END INTERFACE CONTAINS REAL*8 FUNCTION dist( a, b ) TYPE(point), INTENT(IN) :: a, b dist = SQRT( (a%x-b%x)**2 + (a%y-b%y)**2 + (a%z-b%z)**2 ) END FUNCTION dist END MODULE cartesian program over2 USE cartesian IMPLICIT NONE TYPE(point) :: a, b REAL*8 distance a%x=0.0d0 a%y=0.0d0 a%z=0.0d0 b%x=1.0d0 b%y=1.0d0 b%z=1.0d0 ihlee@cetus0:~/lecture_programming$ a.out 1.732050807568877 FORTRAN STOP ihlee@cetus0:~/lecture_programming$ TYPE COORDS_3D REAL :: x, y, z END TYPE COORDS_3D TYPE(COORDS_3D) :: pt1, pt2 distance = a.dist. b write(6,*) distance stop end program over2
MODULE cartesian TYPE point REAL*8 :: x, y END TYPE point INTERFACE ASSIGNMENT( = ) MODULE PROCEDURE max_point END INTERFACE CONTAINS SUBROUTINE max_point( a, pt ) REAL*8, INTENT(OUT) :: a TYPE(point), INTENT(IN) :: pt a = MAX( pt%x, pt%y ) END SUBROUTINE max_point END MODULE cartesian ihlee@cetus0:~/lecture_programming$ a.out 4.200000000000000 FORTRAN STOP ihlee@cetus0:~/lecture_programming$ program over4 USE cartesian IMPLICIT NONE TYPE(point) :: a = point(1.7d0, 4.2d0) REAL*8 :: coord coord = a write(6,*) coord stop end program over4
PROGRAM test INTERFACE REAL*8 FUNCTION ratio(x, y) REAL*8 ::x, y END FUNCTION ratio END INTERFACE ihlee@cetus0:~/lecture_programming$ a.out The ratio is 0.1200000000000000 ihlee@cetus0:~/lecture_programming$ INTEGER :: i=3, j=25 PRINT *,'The ratio is ',ratio(dfloat(i),dfloat(j)) END PROGRAM test REAL*8 FUNCTION ratio(x, y) REAL*8 :: x, y ratio=x/y END FUNCTION ratio
PROGRAM test!interface! REAL*8 FUNCTION ratio(x, y)! REAL*8 ::x, y! END FUNCTION ratio!end INTERFACE ihlee@cetus0:~/lecture_programming$ a.out The ratio is 0.1200000000000000 ihlee@cetus0:~/lecture_programming$ INTEGER :: i=3, j=25 real*8 ratio PRINT *,'The ratio is ',ratio(dfloat(i),dfloat(j)) END PROGRAM test REAL*8 FUNCTION ratio(x, y) REAL*8 :: x, y ratio=x/y END FUNCTION ratio Procedure interface.
SUBROUTINE invert(a, inverse, count) REAL, INTENT(IN) :: a REAL, INTENT(OUT) :: inverse INTEGER, INTENT(INOUT) :: count inverse = 1/a count = count+1 END SUBROUTINE invert PROGRAM test INTERFACE REAL FUNCTION func( x ) REAL, INTENT(IN) ::x END FUNCTION func END INTERFACE... CALL sub1( a, b, func(2) )... END PROGRAM test REAL FUNCTION func( x )! external REAL, INTENT(IN) :: x func = 1/x END FUNCTION func
MODULE strings INTERFACE OPERATOR ( / ) MODULE PROCEDURE num END INTERFACE CONTAINS INTEGER FUNCTION num( s, c ) CHARACTER(len=*), INTENT(IN) :: s CHARACTER, INTENT(IN) :: c num = 0 DO i=1,len( s ) IF( s(i:i)==c ) num=num+1 END DO END FUNCTION num END MODULE strings program ov1 USE strings implicit none character*32, ch1,ch2,ch3 integer ii,jj,kk ihlee@cetus0:~/lecture_programming$ a.out 3 1 1 FORTRAN STOP ihlee@cetus0:~/lecture_programming$ Operator overloading ch1='hello world' ch2='my name is in-ho lee' ch3='who are you' ii = ch1/'l' jj = ch2/'o' kk = ch3/'a' write(6,*) ii write(6,*) jj write(6,*) kk stop end program ov1 http://www.qmw.ac.uk/~cgaa260/building/intr_f90/prg_unts/opoload.htm
MODULE SWAPPER INTERFACE SWAP MODULE PROCEDURE SWAP_R, SWAP_I, SWAP_C END INTERFACE CONTAINS SUBROUTINE SWAP_R(A, B) IMPLICIT NONE REAL, INTENT (INOUT) :: A, B REAL :: TEMP TEMP = A ; A = B ; B = TEMP END SUBROUTINE SWAP_R ihlee@cetus0:~/lecture_programming$ a.out 1 2 100 200 7.100000 10.90000 11.10000 17.00000 ab1" 2 1 200 100 10.90000 7.100000 17.00000 11.10000 ba"1 ihlee@cetus0:~/lecture_programming$ SUBROUTINE SWAP_I(A, B) IMPLICIT NONE INTEGER, INTENT (INOUT) :: A, B INTEGER :: TEMP TEMP = A ; A = B ; B = TEMP END SUBROUTINE SWAP_I SUBROUTINE SWAP_C(A, B) IMPLICIT NONE CHARACTER, INTENT (INOUT) :: A, B CHARACTER :: TEMP TEMP = A ; A = B ; B = TEMP END SUBROUTINE SWAP_C END MODULE SWAPPER PROGRAM SWAP_MAIN USE SWAPPER IMPLICIT NONE INTEGER :: I, J, K, L REAL :: A, B, X, Y CHARACTER :: C, D, E, F I = 1 ; J = 2 ; K = 100 ; L = 200 A = 7.1 ; B = 10.9 ; X = 11.1; Y = 17.0 C = 'a' ; d = 'b' ; E = '1' ; F = '"' WRITE (*,*) I, J, K, L, A, B, X, Y, C, D, E, F CALL SWAP (I, J) ; CALL SWAP (K, L) CALL SWAP (A, B) ; CALL SWAP (X, Y) CALL SWAP (C, D) ; CALL SWAP (E, F) WRITE (*,*) I, J, K, L, A, B, X, Y, C, D, E, F END http://www-teaching.physics.ox.ac.uk/unix+prog/hargrove/tutorial_90/13_extra.html
Derived Type TYPE COORDS_3D REAL :: x, y, z END TYPE COORDS_3D TYPE(COORDS_3D) :: pt1, pt2 TYPE SPHERE TYPE (COORDS_3D) :: centre REAL :: radius END TYPE SPHERE TYPE (SPHERE) :: bubble, ball pt1%x = 1.0 bubble%radius = 3.0 bubble%centre%x = 1.0 pt1 = COORDS_3D(1.,2.,3.) bubble%centre = COORDS_3D(1.,2.,3.) bubble = SPHERE(bubble%centre,10.) bubble = SPHERE(COORDS_3D(1.,2.,3.),10.)
Multiple Precision Computation. 2**200 = 1606938044258990275541962092341162602522202993782792835301376 (prime number)? 5468317884572019103692012212053793153845065543480825746529998049913561 The same example using derived types:... USE FMZM TYPE ( FM ) A,B,C,D... D = SQRT( A**2 + B**2 + C**2 )... http://incredible.egloos.com/3160240
References http://incredible.egloos.com/2950367 pdf.. 90. http://www.fortran.com/fortran/ Pointer to everything Fortran http://www.arc.unm.edu/~jabed/list.html Pointer to a list of tutorials http://dynaweb.sdsc.edu:8080/library/t90_c90_y90_cray_fp?dweb_collection=dweb_collections Cray's Manual http://www.software.ibm.com/ad/fortran/xlfortran/optim.html Short optimization guide from IBM http://www.digital.com/fortran/dvf.html DEC's Digital Visual Fortran Page ftp://mycroft.plk.af.mil/pub/fortran_90/tutorial A tutorial by Zane Dodson http://www.nova.edu/ocean/psplot.html Postscript plotting library http://meteora.ucsd.edu/~pierce/fxdr_home_page.html Subroutines to do unformatted I/O across platforms, provided by David Pierce at UCSD http://www.nsc.liu.se/~boein/f77to90/a5.html A good reference for intrinsic functions http://www.nsc.liu.se/~boein/f77to90/ Fortran 90 for the Fortran 77 Programmer http://www.tat.physik.uni-tuebingen.de/~kley/lehre/ftn77/tutorial/index.html http://www.sdsc.edu/~tkaiser/f90.html
(?)
Very-high-level language (shell) (very-high-level language). (.). (!). (; ; object-oriented)..... No compile Monty Python's Flying Circus. No type declarations No memory allocation control High level data structure : (Guido van Rossum) http://en.wikibooks.org/wiki/python_programming http://incredible.egloos.com/2950369 http://en.wikipedia.org/wiki/python_(programming_language) http://pentangle.net/python/handbook/handbook.html http://www.pasteur.fr/recherche/unites/sis/formation/python/
. (dynamic typing),., (),,. (glue language)..,. (dynamic typing). (.).,,. ("Battery Included"),.. (MS windows,,,.,.)! : http://www.awaretek.com/tutorials.html#begin http://ko.wikipedia.org/wiki/%ed%8c%8c%ec%9d%b4%ec%8d%ac
Python Control-d exit(). ihlee@cetus0:~/lecture_programming$ python Python 2.3.5 (#2, Oct 16 2006, 19:59:54) [GCC 3.3.5 (Debian 1:3.3.5-13)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> print "Hello world" ihlee@cetus0:~/lecture_programming$ python Python 2.3.5 (#2, Oct 16 2006, 19:59:54) [GCC 3.3.5 (Debian 1:3.3.5-13)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> a=1. ; b=2.d0 File "<stdin>", line 1 a=1. ; b=2.d0 ^ SyntaxError: invalid syntax >>> a=1 ; b=2. >>> c=a+b >>> print c 3.0 >>> bc l. quit. ihlee@cetus0:~$ bc -l bc 1.06 Copyright 1991-1994, 1997, 1998, 2000 Free Software Foundation, Inc. This is free software with ABSOLUTELY NO WARRANTY. For details type `warranty'. 12.*12 144 quit ihlee@cetus0:~$
* (indentation) python.! *# comment. ``. * \.) "". * (,,,. *. #!/usr/bin/env python. * "def" ":".. if for. *for i in range(100): i 0 0,1,2,3,4,...99 for loop 100. range range(1,101) 1 100 i. range(10) ( ). [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] *FORTRAN 90.. * return.,... return. return a,b,c,d. a, b 1, c 2 ( ) d.. abcd() aa,bb,cc,dd=abcd().
#!/usr/bin/env python #from Numeric import * #from string import atof,atoi,split #from posix import getpid,rename,remove,mkdir,rmdir,system #from shutil import copyfile #import os,sys,time,random version ="1.0.0" counter=0 acc = 0 for n in xrange(1,10000000000): num1 = str(n).count('5') if num1 > 0 : acc = acc+num1 if acc == n: counter=counter+1 print "counter = ",counter,n if counter == 2 : break,. chmod 755 abc.py
system( )... ` `;? CALL SYSTEM('application.py'). call system. ( ).... : / == call system('python_driver.py') ==. : / ==
sad.x 2000.,.. for i in range(2000): # 2000 gen_sad_input() # sad.i system('sad.x >TMPFILE') # sad.x copyfile('dim.dat', ascnum('dimer.',i)) copyfile('mid.xyz', ascnum_xyz('midpt.',i)) # system('rm dim.dat') # system('rm TMPFILE') os.path.exists >>> from os import * >>> system('ls') abc.py HFB PBS_script_examples tr.py action_protein_draco JOB_ID_PRINTED swkim_mpi VASP_example admd_cluster matplot_test test v_tubes dssp onion test1 ykm 0 >>> system(), system('qsub work.0003/pbs_script_file') pbs. j=3 ; f_name=ascnum('work.',j)+'/pbs_script_file' ; system('qsub '+f_name+'') import glob,time import os.path file_list=glob.glob('*') for f_name in file_list: if os.path.isfile(f_name): print f_name,' is a regular file' print os.path.getsize(f_name) print os.path.getatime(f_name) print time.ctime(os.path.getatime(f_name)) elif os.path.isdir(f_name): print f_name,' is a directory' elif os.path.islink(f_name): print f_name,' is a symbolic link'
#!/usr/bin/env python faren_old = float(raw_input("enter yesterday's temperature ")) faren_new = float(raw_input("enter today's temperature ")) if faren_new > 90 and faren_old < 90: print "First day of a heatwave" elif faren_new > 90: print "Subsequent day of a heatwave" elif faren_old < 90: print "Not the hot season" else: print "Bye, bye heatwave" print "Program completed" >>> raw_input() 3 '3' >>> input() 4 4 ihlee@cetus0:~/lecture_programming$ py1.py Enter yesterday's temperature 99 Enter today's temperature 98. Subsequent day of a heatwave Program completed ihlee@cetus0:~/lecture_programming$
phi_trj=[[0.]*(iidd) for j in xrange(np+1)] psi_trj=[[0.]*(iidd) for j in xrange(np+1)] for ifile in range(0,np+1): tmp_file=asc2num('rama',ifile)+'.txt' fp_rama=open(tmp_file,'r') jth=0 while 1: line=fp_rama.readline() if not line : break else: token=line.split() phi_trj[ifile][jth]=float(token[0]) psi_trj[ifile][jth]=float(token[1]) jth=jth+1 fp_rama.close() http://incredible.egloos.com/3588647
def wwait(file_name): while 1 : time.sleep(60.) if os.path.exists(file_name) : remove(file_name) break return # 60 # def mkdir_series(n_dir): for j in range(n_dir): char=ascnum('test_dir.',j) # test_dir.0001 n_dir mkdir(char) return touch abcd
#!/usr/bin/env python from Numeric import * from string import atof,atoi,split from posix import getpid,rename,remove,mkdir,rmdir,system from shutil import copyfile import os,sys,time,random version ="1.0.0" if name == " main ": # # *.bin # convert.x #. convert.x "convert.x test.bin test". # # yeskim98, March/26/2003 # system('ls *.bin >namelist') # f=open('namelist','r') for line in f: name=line.split()[0] name1=name[:-4] # 4 system('./convert.x '+name+' '+name1+'') f.close()
#!/usr/bin/env python from Numeric import * # from string import atof,atoi,split from posix import getpid,rename,remove,mkdir,rmdir,system from shutil import copyfile import os,sys,time,random version ="1.0.0" def prtjobid(): # written by In-Ho Lee, KRISS, Jan. 12, 2002. f=open('job_id_printed','w') f.write('%15d \n' % os.getpid()) # st=time.asctime() # f.write('%23s \n' % st) f.close() return def num2char(i): # ic=str(i) ic='0' if i == 1 : ic='1' if i == 2 : ic='2' if i == 3 : ic='3' if i == 4 : ic='4' if i == 5 : ic='5' if i == 6 : ic='6' if i == 7 : ic='7' if i == 8 : ic='8' if i == 9 : ic='9' return ic def ascnum(ichar,number): # written by In-Ho Lee, KRISS, Jan. 12, 2002. i=int(float(number)/1000.) j=int((number-i*1000.)/100.) k=int((number-j*100.-i*1000.)/10.) l=int(number-k*10.-j*100.-i*1000.) ic=num2char(i) jc=num2char(j) kc=num2char(k) lc=num2char(l) fname=ichar+ic+jc+kc+lc # abc.0001 return fname ihlee@cetus0:~/lecture_programming$ num.py abc0000.xyz abc0001.xyz abc0002.xyz abc0003.xyz abc0004.xyz abc0005.xyz abc0006.xyz abc0007.xyz abc0008.xyz abc0009.xyz ihlee@cetus0:~/lecture_programming$ >>> a=1 >>> b=str(a).zfill(4) >>> print b 0001 >>> c=str(a).zfill(5) >>> print c 00001 >>> c=str(c) >>> print c 00001 >>> >>> from string import * >>> a=1 >>> b=str(a) >>> print b 1 >>> a=1.e0 >>> c=str(222) >>> a=complex(a) >>> print c >>> print a 222 (1+0j) >>> >>> def ascnum_xyz(ichar,number): # written by In-Ho Lee, KRISS, Jan. 12, 2002. fname=ascnum(ichar,number)+'.xyz' # abc.0001.xyz return fname prtjobid() ichar='abc' for i in range(10): print ascnum_xyz(ichar,i) >>> print atof(c) 222.0 >>> print atoi(c) 222 >>> print float(c) 222.0 >>> print int(c) 222
Prime number? #!/usr/bin/env python for i in range(2,10): for x in range(2,i): if i % x == 0 : print i, 'equals', x, '*', i/x break else: print i, 'is a prime number' ihlee@cetus0:~/lecture_programming$ prime.py 3 is a prime number 4 equals 2 * 2 5 is a prime number 5 is a prime number 5 is a prime number 6 equals 2 * 3 7 is a prime number 7 is a prime number 7 is a prime number 7 is a prime number 7 is a prime number 8 equals 2 * 4 9 is a prime number 9 equals 3 * 3 ihlee@cetus0:~/lecture_programming$