Transcript AME 150 L - Engineering Class Home Pages

AME 150 L
Fortran Statements
Functions & Simple IO
7 - 1/26/2000
AME 150 L
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
7 - 1/26/2000
AME 150 L
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
7 - 1/26/2000
AME 150 L
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
7 - 1/26/2000
AME 150 L
Role of Declaration
Declaration statement tells compile what type
of arithmetic to use with variable
Operand
Other Operand
Arithmetic
INTEGER
INTEGER
INTEGER
INTEGER
REAL
REAL
REAL
REAL
REAL
7 - 1/26/2000
AME 150 L
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 - 1/26/2000
AME 150 L
Mixed Mode Examples
INTEGER :: I, J, K
REAL X, Y, Z
I+J ?
Integer
X+Y ?
Real
I+X ?
Real
7 - 1/26/2000
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)
AME 150 L
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?)
7 - 1/26/2000
AME 150 L
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 + -
7 - 1/26/2000
AME 150 L
Arithmetic Types
• Arithmetic operation affect ONLY
arithmetic variables
– INTEGER
– REAL
– COMPLEX
– Multiple precision (machine dependent)
– Special TYPE based on arithmentic
7 - 1/26/2000
AME 150 L
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
7 - 1/26/2000
AME 150 L
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.
7 - 1/26/2000
AME 150 L
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
7 - 1/26/2000
AME 150 L
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
Octal
HexaDecimal
7 - 1/26/2000
0.33333334
7652525253
3EAAAAAB
0.25000000
7640000000
3E800000
AME 150 L
0.20000000
7623146315
3E4CCCCD
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 7 - 1/26/2000
AME 150 L
Examples
• Integer
13 = 23+22+20= 11012=158=d16
• Real
13.0=23*(20+2-1+2-3)exponent=8=10002
and 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
7 - 1/26/2000
AME 150 L
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
7 - 1/26/2000
AME 150 L
Sample Calculation
• Given: 3 coefficients a,b,c of quadratic
equation: ax2+bx+c=0
• Find 2 roots of the equation x1, x2 where
2

b

b  4ac
x1,2 
2a
• More in B: p87-91, D: p156-161
7 - 1/26/2000
AME 150 L
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
7 - 1/26/2000
AME 150 L
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
7 - 1/26/2000
AME 150 L
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
7 - 1/26/2000
AME 150 L
Trial 2
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
7 - 1/26/2000
AME 150 L
!These lines will
!be replaced
Creating T2.f90
PROGRAM T2
REAL :: a,b,c,x1,x2,disc
READ(*,*)a,b,c
disc = b**2 - 4. * a * c
IF (disc >= 0 ) THEN
x1 = ( -b + sqrt(disc) )/(2.*a)
x2 = ( -b - sqrt(disc) )/(2.*a)
ELSE
!Do something
END IF
WRITE(*,*)"x1,x2=",x1,x2
END PROGRAM T2
7 - 1/26/2000
!declare a new variable
!calculate this 1 time
!Executed when (…) is .TRUE.
! The ELSE clause is executed
!when (…) is .FALSE.
AME 150 L
Complex Conjugate Roots
• When disc < 0, then calculate i*sqrt(-disc)
– where i2 = -1
• Complex Conjugate Roots
Real Part:
 Imaginary part
7 - 1/26/2000
xreal = -b/(2.*a)
ximag= SQRT(-disc)
AME 150 L
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
7 - 1/26/2000
AME 150 L
Special Case a=0
• Test first for a /= 0, and do "t2.f90"
calculations if .TRUE.
• IF (a==0) [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)
7 - 1/26/2000
AME 150 L
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
7 - 1/26/2000
AME 150 L
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
7 - 1/26/2000
AME 150 L
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
7 - 1/26/2000
AME 150 L
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
7 - 1/26/2000
AME 150 L
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
7 - 1/26/2000
AME 150 L
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?
7 - 1/26/2000
AME 150 L
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
7 - 1/26/2000
AME 150 L