1 / 33

# AME 150 L

AME 150 L. Fortran Statements Functions &amp; 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

## AME 150 L

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