Fortran Statements

1
AME 150 L

• Fortran Statements
• Functions Simple IO

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

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

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

5
Role of Declaration

• Declaration statement tells compile what type of
  arithmetic to use with variable

6
Mixed Mode Arithmetic

• Mixed Mode results for every expression takes
  form of most complex
• Avoid if possible (confusion can result)
• Never mix CHARACTER or LOGICAL with REAL or
  INTEGER
• Problems and confusion occurs more with constants
  than with variables

7
Mixed Mode Examples

• INTEGER I, J, K
• REAL X, Y, Z
• IJ ?
• Integer
• XY ?
• Real
• IX ?
• 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)

8
Exponentiation

• Real Integer ? Repeated Multiplication
  (-7.53)
• Real -Integer ? 1./(Real Integer )
• Integer Integer ? Repeated Multiplication
  (57)
• Integer -Integer ? 0 (Why?)
• I-J 1/(IJ) ? 0 But K /IJ not necessarily
  0
• Integer Real ? Never Allowed (Why?)
• X Y ? Exponential (Y log(X))
• X MUST BE Positive (Why?)

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 -

10
Arithmetic Types

• Arithmetic operation affect ONLY arithmetic
  variables
• INTEGER
• REAL
• COMPLEX
• Multiple precision (machine dependent)
• Special TYPE based on arithmentic

11
Other Types

• LOGICAL
• Special Operators (AND, OR, NOT, etc)
• Can create by relational operations
• X gt Y (read X greater than Y) is ONLY .True. or
  .False.
• Relational Operators
• Equal to / Not equal to
• gt Greater than lt Less than or equal to
• lt Less than gt Greater than or equal to

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.

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

14
PROGRAM report implicit none REAL third,
fourth,fifth third1./3. fourth1./4. fifth1./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
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 1fraction 2exponent
• fraction is less than 1, but number is 1.bbbb
• both fraction exponent can be or -

16
Examples

• Integer
• 13 232220 11012158d16
• Real
• 13.023(202-12-3)?exponent810002and the
  number 1.5.1251.625101.1012
• Negative numbers -- use 2's complement
• Change all 1s to 0s and all 0s to 1's (ones
  complement) and then add 1

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

18
Sample Calculation

• Given 3 coefficients a,b,c of quadratic
  equation ax2bxc0
• Find 2 roots of the equation x1, x2 where
• More in B p87-91, D p156-161

19
Trial 1

• PROGRAM T1
• REAL a,b,c,x1,x2
• READ(,)a,b,c
• x1(-bsqrt(b2-4.ac))/(2.a)
• x2(-b-sqrt(b2-4.ac))/(2.a)
• WRITE(,)"x1,x2",x1,x2
• END PROGRAM T1

20
Problems with T1

• 1) a0 -- divide by 0
• must test for a0, and if so, the equation is
• bxc0, so that only one solution exists x-c/b
  (as long as b ?0 -- not important since c0 too)
• 2) a ?0, but b2-4ac lt0, so the two roots are a
  complex conjugate pair
• 3) a ?0, and b2-4ac 0, so there are two equal
  roots

21
First Repair

• Fix the b2-4ac lt0 problem first
• b2-4ac (and its square root) is computed twice
• Let's test for negative first
• Do an intermediate calculation
• disc b2-4.ac
• and make sure disc gt0 before taking SQRT
• Need something like "if disc is less than 0,
  then"
• IF (disc lt 0) THEN

22
Trial 2

• PROGRAM T1
• REAL a,b,c,x1,x2
• READ(,)a,b,c
• x1(-bsqrt(b2-4.ac))/(2.a) !These lines
  will
• x2(-b-sqrt(b2-4.ac))/(2.a) !be replaced
• WRITE(,)"x1,x2",x1,x2
• END PROGRAM T1

23
Creating T2.f90

• PROGRAM T2
• REAL a,b,c,x1,x2,disc !declare a new variable
• READ(,)a,b,c
• disc b2 - 4. a c !calculate this 1 time
• IF (disc gt 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

24
Complex Conjugate Roots

• When disc lt 0, then calculate isqrt(-disc)
• where i2 -1
• Complex Conjugate Roots
• Real Part xreal -b/(2.a)
• ? Imaginary part ximag SQRT(-disc)

25
Finishing T2.f90

• REAL a,b,c,x1,x2,disc,xreal,ximag !declare a
  new variables
• disc b2 - 4. a c !calculate this 1 time
• IF (disc gt 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

26
Special Case a0

• Test first for a / 0, and do "t2.f90"
  calculations if .TRUE.
• IF (a0) i.e. (a/ 0) .FALSE. Use ELSE
• Write a message
• Calculate one root (answer -c/b)
• Save these as t3.f90
• t4.f90 is cleaned up (and no extra calculations)

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,c0.333333,.666667,.333333 had roots x1,
  x2 -0.998264432, -1.00173867
• ARITHMETIC ROUND-OFF

28
No Logical Close-to

• disc0 is hard to make happen
• REAL number calculation tries its best to round
  off, but it is hard
• Sample error - counting with REAL

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 xmin0 xmax100. deltax(xmax-xmin)/nste
ps Count0 xxmin DO xxdeltax
CountCount1 IF (x xmax)EXIT END
DO WRITE(,)x, Count Stop END PROGRAM xcount
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.3333333 0.999999 -- close to, but not 1
• The relation real1 real2 is a problem

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_epsilon1.19209E-07

32
Using Machine Epsilon

• Differences in real numbers
• x1-x2 lt eps
• Need the absolute value of the difference
• ABS(x1-x2) lt eps
• What happens if x1,x2 is very large (small)
• ABS(x1-x2) / ABS(x1x2) lt eps
• Is this ok?

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) lt epsABS(x1x2)
• Almost, but certainly, both x1 and x2 0 should
  pass this test
• ABS(x1-x2) lt epsABS(x1x2)
• is .TRUE. when x1 is approximately x2