1 / 33

AME 150 L

AME 150 L. Fortran Statements Functions & Simple IO. Elementary Input/Output. Input: READ ( unit, format ) list of variables Simple form: READ (*,*) list Read from Standard Input (keyboard unless redirected) Output: WRITE ( unit, format ) list of variables

zubin
Download Presentation

AME 150 L

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. AME 150 L Fortran Statements Functions & Simple IO AME 150 L

  2. Elementary Input/Output • Input:READ (unit, format) list of variables • Simple form: READ (*,*)list • Read from Standard Input (keyboard unless redirected) • Output: WRITE(unit, format) list of variables • Simple form: WRITE (*,*)list • Display on Standard Output (display unless redirected) • Reading &/or writing is done upon execution AME 150 L

  3. Model Program PROGRAM model1 IMPLICIT NONE ! Declarations of input data & results READ ( *,*) input data results = calculations on input data WRITE (*,*) results STOP END PROGRAM model1 AME 150 L

  4. Declarations • IMPLICIT NONE ALL variables must be declared • Declared means that the variable name must appear in an INTEGER or REAL or other TYPE statement, i.e. REAL :: velocity, a1, x, y, z, Mach_Number INTEGER :: index, I, j, k, Counter AME 150 L

  5. Role of Declaration Declaration statement tells compile what type of arithmetic to use with variable AME 150 L

  6. Mixed Mode Arithmetic • Mixed Mode results for every expression takes form of most complex • Avoid if possible (confusion can result) • Never mixCHARACTER or LOGICAL with REAL or INTEGER • Problems and confusion occurs more with constants than with variables AME 150 L

  7. INTEGER :: I, J, K REAL X, Y, Z I+J ? Integer X+Y ? Real I+X ? Real I / J ? Integer 5 / 2 ? Integer = 2 (not 2.5) X + I / K ? Real + Integer = Real 2.7 + 5 / 2 ? 2.7 + 2 = 4.7 (not 5.2) Mixed Mode Examples AME 150 L

  8. Exponentiation Real ** Integer  Repeated Multiplication (-7.5**3) Real ** -Integer  1./(Real ** Integer ) Integer ** Integer  Repeated Multiplication (5**7) Integer ** -Integer  0 (Why?) I**-J = 1/(I**J)  0 But K /I**J not necessarily 0 Integer ** Real  Never Allowed (Why?) X ** Y  Exponential (Y * log(X)) X MUST BE Positive (Why?) AME 150 L

  9. Reminder - Hierarchy of Operations • First, evaluate all (…) including functions • Start evaluations without (..) with • Unary minus - • Exponentiation ** • Multiply & Divide (left to right generally) * / • Add and Subtract + - AME 150 L

  10. Arithmetic Types • Arithmetic operation affect ONLY arithmetic variables • INTEGER • REAL • COMPLEX • Multiple precision (machine dependent) • Special TYPE based on arithmentic AME 150 L

  11. Other Types LOGICAL • Special Operators (AND, OR, NOT, etc) • Can create by relational operations X > Y(read X greater than Y) is ONLY .True. or .False. Relational Operators == Equal to /= Not equal to > Greater than <= Less than or equal to < Less than >= Greater than or equal to AME 150 L

  12. Relational Operators • Can mix mode, but must be careful • Can use them on CHARACTER data type • Cannot mix Character with Arithemetic • Always .true. or .false. Never .maybe. • Relational operators are lower in the hierarchy than all arithmetic • There is no .approximately. or .close to. AME 150 L

  13. The Real Curse • All computers use binary • We use decimal • We know 1./3. = 0.333333333… • and 1./5. = 0.2 (exact) • In Binary, both repeat indefinitely AME 150 L

  14. PROGRAM report implicit none REAL :: third, fourth,fifth third=1./3. fourth=1./4. fifth=1./5. WRITE(*,*)third, fourth,fifth WRITE(*,"(3o20)" ) third, fourth,fifth stop END PROGRAM report Output: Decimal 0.33333334 0.25000000 0.20000000 Octal 7652525253 7640000000 7623146315 HexaDecimal 3EAAAAAB 3E800000 3E4CCCCD AME 150 L

  15. Storage of REAL Variables • Engineering or Science notation, 1.23 105 • 1.23 has significant figures (in decimal) • 5 is the exponent (power of 10) • Normalized Floating Point 0.123 106 • fraction is between 1 and 1/10th • Computer data is binary 1+fraction 2exponent • fraction is less than 1, but number is 1.bbbb • both fraction & exponent can be + or - AME 150 L

  16. Examples • Integer 13 = 23+22+20= 11012=158=d16 • Real 13.0=23*(20+2-1+2-3)exponent=8=10002and the number =1+.5+.125=1.62510=1.1012 • Negative numbers -- use 2's complement • Change all 1s to 0s and all 0s to 1's (ones complement) and then add 1 AME 150 L

  17. IEEE Floating Point Numbers • I found a great web site -- it's from a Florida State University Numerical Analysis Class http://www.scri.fsu.edu/~jac/MAD3401/Backgrnd/ieee.html • I have placed this link in the "useful stuff" page on the AME150 web page AME 150 L

  18. Sample Calculation • Given: 3 coefficients a,b,c of quadratic equation: ax2+bx+c=0 • Find 2 roots of the equation x1, x2 where • More in B: p87-91, D: p156-161 AME 150 L

  19. Trial 1 PROGRAM T1 REAL :: a,b,c,x1,x2 READ(*,*)a,b,c x1=(-b+sqrt(b**2-4.*a*c))/(2.*a) x2=(-b-sqrt(b**2-4.*a*c))/(2.*a) WRITE(*,*)"x1,x2=",x1,x2 END PROGRAM T1 AME 150 L

  20. Problems with T1 1) a=0 -- divide by 0 must test for a=0, and if so, the equation is bx+c=0, so that only one solution exists x=-c/b (as long as b 0 -- not important since c=0 too) 2) a 0, but b2-4ac <0, so the two roots are a complex conjugate pair 3) a 0, and b2-4ac =0, so there are two equal roots AME 150 L

  21. First Repair • Fix the b2-4ac <0 problem first b2-4ac (and its square root) is computed twice Let's test for negative first Do an intermediate calculation disc = b**2-4.*a*c and make sure disc >0 before taking SQRT Need something like "if disc is less than 0, then" IF (disc < 0) THEN AME 150 L

  22. Trial 2 PROGRAM T1 REAL :: a,b,c,x1,x2 READ(*,*)a,b,c x1=(-b+sqrt(b**2-4.*a*c))/(2.*a)!These lines will x2=(-b-sqrt(b**2-4.*a*c))/(2.*a) !be replaced WRITE(*,*)"x1,x2=",x1,x2 END PROGRAM T1 AME 150 L

  23. PROGRAM T2 REAL :: a,b,c,x1,x2,disc !declare a new variable READ(*,*)a,b,c disc = b**2 - 4. * a * c !calculate this 1 time IF (disc >= 0 ) THEN x1 = ( -b + sqrt(disc) )/(2.*a) !Executed when (…) is .TRUE. x2 = ( -b - sqrt(disc) )/(2.*a) ELSE ! The ELSE clause is executed !Do something !when (…) is .FALSE. END IF WRITE(*,*)"x1,x2=",x1,x2 END PROGRAM T2 Creating T2.f90 AME 150 L

  24. Complex Conjugate Roots • When disc < 0, then calculate i*sqrt(-disc) • where i2= -1 • Complex Conjugate Roots Real Part: xreal = -b/(2.*a)  Imaginary part ximag= SQRT(-disc) AME 150 L

  25. Finishing T2.f90 REAL :: a,b,c,x1,x2,disc,xreal,ximag !declare a new variables … disc = b**2 - 4. * a * c !calculate this 1 time IF (disc >= 0 ) THEN … WRITE(*,*)"x1,x2=",x1,x2 !move the Real root write ELSE ! The ELSE clause is executed xreal = -b/(2.*a) !calculate real part ximag = sqrt(-disc)/(2.*a) !and imaginary part WRITE(*,*)"Complex Conjugate roots",xreal,"+-I*",ximag END IF END PROGRAM T2 AME 150 L

  26. Special Case a=0 • Test first for a /= 0, and do "t2.f90" calculations if .TRUE. • IF (a==0)[i.e. (a/= 0) .FALSE.UseELSE] Write a message Calculate one root (answer = -c/b) • Save these as t3.f90 • t4.f90 is cleaned up (and no extra calculations) AME 150 L

  27. Weird Problem • Case a,b,c = 1,2,1 has roots -1,-1 • Case a,b,c = .33333,.66666,.33333 also has roots -1,-1 • BUT a,b,c=.66667, 1.33333,.66667 caused an arithmetic exception and dumped core • AND a,b,c=0.333333,.666667,.333333 had roots x1, x2 = -0.998264432, -1.00173867 ARITHMETIC ROUND-OFF AME 150 L

  28. No Logical Close-to • disc==0 is hard to make happen • REAL number calculation tries its best to round off, but it is hard • Sample error - counting with REAL AME 150 L

  29. PROGRAM xcount IMPLICIT none REAL :: x,xmin,xmax,deltax,nsteps INTEGER :: Count WRITE(*,*)"How many steps do you want?" READ(*,*)nsteps !Initialize the variables xmin=0 xmax=100. deltax=(xmax-xmin)/nsteps Count=0 x=xmin DO x=x+deltax Count=Count+1 IF (x == xmax)EXIT END DO WRITE(*,*)x, Count Stop END PROGRAM xcount AME 150 L

  30. What is happening? • Sometimes the program works, sometimes it doesn't • All failures occur for cases where (xmax-xmin)/nsteps is a repeating binary number (like 0.333333 is for 1/3.) 0.333333*3 =0.999999 -- close to, but not == 1 • The relation real1 == real2 is a problem AME 150 L

  31. Machine Epsilon • Epsilon (): Math. a small, positive number • Machine epsilon refers to the least significant bit -- often it is 2-precision in REAL numbers on your computer • Example: If there is 23 bits in a REAL fraction, machine_epsilon=1.19209E-07 AME 150 L

  32. Using Machine Epsilon • Differences in real numbers x1-x2 < eps • Need the absolute value of the difference ABS(x1-x2) < eps • What happens if x1,x2 is very large (small) ABS(x1-x2) / ABS(x1+x2) < eps • Is this ok? AME 150 L

  33. Machine Epsilon (continued) • Problem is, there is a divide (slow), and if both x1 & x2 are 0, a divide by 0 error ABS(x1-x2) < eps*ABS(x1+x2) • Almost, but certainly, both x1 and x2 =0 should pass this test ABS(x1-x2) <= eps*ABS(x1+x2) is .TRUE. when x1 is approximately= x2 AME 150 L

More Related