= .TRUE.

Chapter 3 e.g. Mixed operations: (a+b) < (c+d) .OR. 3>=4

Precedence order

: 1. arithmetic operators (*, +, etc, evaluated left-to right) ; 2. relational operators ( < , = =, etc, evaluated from left to right) 3. .NOT. 4. .AND. (evaluated from left to right) 5. .OR. (evaluated from left to right) 6. .EQV. and .NEQV. (evaluated from left to right)

Chapter 3 Note: combinational operators are illegal to use with real/integer/character e.g. 4 .AND. 3 (compile-error)

3.3.5 Logical values for Input/Output statements

When entering data, using READ(*,*), if the first entered character is ‘T’ then .TRUE. is read in

;

similarly for .FALSE.

3.3.6 Significance of Logical variables/expressions

They are mostly used with IF statements.

Chapter 3

3.4 Control constructs: Branches Branche

s in Fortran allow us to execute (depending on the value true/false of logical values/expressions) specific parts of the program and skip other parts. The block IF construct is: IF (logical_expression )

THEN ELSE IF

(

logical

_

expression

)

THEN

IF(.NOT. (a < b. AND. b < c)) ! Example of block IF x = a + b ! e.g. a=3,b=4, c=5 ELSE IF (a < b. AND. b < c) x = c + b ENDIF

Chapter 3 ,

3.4.2 The ELSE and ELSE IF 3.4.2 The ELSE and ELSE IF

:

IF

(criterion_1)

THEN

action_1

ELSE IF (criterion_2) THEN

action_2

ELSE IF (criterion_3) THEN

action_3

ELSE END IF

action_4

A minimal block

IF: IF

(

criterion

)

THEN

action

END IF

Chapter 3

Example of Fortran program:

The quadratic equation (

PROGRAM roots, Fig-3-10

). Examples of inputs to test the program: x^2+5x+6 =0, roots: -2, -3 x^2+4x+4 =0, roots: -2 (double) x^2+2x+5 =0, roots: -1+2i, -1-2i (complex)

Chapter 3

Example of Fortran program:

Evaluating a function of two variables (

PROGRAM funxy, Fig-3-12

). Test the program: f(2,3) = 5 f(-2,3) = 7 f(2,-3) = 11 f(-2,-3) = 13

Chapter 3

3.4.4 The named block IF constructs

:

[name:] IF

(

logical expression

)

THEN

statement1

ELSE IF

(

logical expression

)

THEN

statement1

ELSE

statement1

END IF

Chapter 3

3.4.5 Notes concerning use of block IF constructs

: Block IF may be

nested: outer: IF

(x>0.)

THEN

statement1

inner: IF

(y<0.)

THEN

statement1

END IF inner

statement1

END IF outer

Chapter 3

3.4.6

The logical IF statement: e.g. IF(x<0.) y=-1

Note

: Better to avoid it and use Block IF, instead.

3.4.7 The SELECT CASE construct

deals with many alternatives are mutually exclusive. The case selector must be : integer, character, or logical expression. But ‘real’ expressions are not allowed.

Chapter 3

[name:]SELECT CASE

(

case expression

)

CASE

(

case selector

) !mutually exclusive case

block of Fortran statements

CASE

(

case selector

) !mutually exclusive case

block of Fortran statements

.

.

END SELECT name

Chapter 3

Example of Fortran program:

Use SELECT CASE to select the day of the week (

PROGRAM day_of_week

, Fig-3-14). Test: e.g. 1 Day= Sunday

Chapter 3

Example of Fortran program:

in a SELECT CASE construct Use characters (

PROGRAM weekday_weekend

, Fig. 3-15). Test: e.g. Tuesday Day Type= weekday

Chapter 3

Example

: Assigning letter grades 90< grade A 80< grade <=89 B 70< grade <=79 C 61< grade <=69 D 0< grade <=60 F

Chapter 3

(a)

Solution without

nested

(start from highest grade and create mutual exclusive logical exressions)

IF

(grade > 90)

THEN

WRITE(*,*) ‘The grade is A‘

ELSE IF

(grade > 80)

THEN

WRITE(*,*) ‘The grade is B‘

ELSE IF

(grade > 70 )

THEN

WRITE(*,*) ‘The grade is C ‘

ELSE IF

(grade > 61 )

THEN

WRITE(*,*) ‘The grade is D‘

ELSE

WRITE(*,*) ‘The grade is F‘

END IF

Chapter 3

(b)

Solution using nested IF (in textbook p. 111), to be avoided because it may lead to errors.

(c)

Solution using SELECT CASE left as an exercise (not to be handed in).

Chapter 3

3.5

Debugging Fortran programs We can

insert

WRITE(*,*) statements into the program to

print out important var’s at key points

program. e.g. Verify correctness of block IF: WRITE(*,*) ‘At if1: var1 =‘, var1 of the if1: IF(sqrt(var1) >1.) THEN WRITE(*,*) ‘At if1: sqrt(var1) >1. ‘ ELSE IF(sqrt(var1) <1.) THEN WRITE(*,*) ‘At if1: sqrt(var1) <1. ‘ ELSE WRITE(*,*) ‘At if1: sqrt(var1) == 1. ‘ END IF if1

Chapter 3 When we locate in which part of the program the error occurs. We check the statements to discover it. Examples of common errors with IF statements: 1.

Check the relational operators e.g. maybe >= is needed instead of > , etc.

2. When we test for equality of real var’s/const’s (unlike integers), we can not use == .

We must use e.g. IF(ABS(x-10.) < 0.000001) THEN IF(x==10.), for x real, may be false because in some computer x=0.9999999 (which is 10.0).

Chapter 3

Recommended Problems (not to be handed in)

• Page 95 , Quiz 3-1 (answers are on the back): 1(a,b,d,f,h), 2 • Page 117 , Quiz 3-2 (answers are on the back): 1,2,4,6,7 Page 123-124, Exercises: 3-1, 3-2, 3-4, 3-6, 3-7

Chapter 4 LOOPS AND CHARACTER MANIPULATION Objectives

• Know how to create and use while loops.

• Know how to create and use counting loops.

• Know when to use while loops, and when to use counting loops. • Know the purpose of the CONTINUE and EXIT statements, and how to use them. • Understand loop names, and why they are used. • Learn about character assignments and character operators. • Learn about substrings and string manipulations.

While Loop

The While Loop

DO ...

IF (

logical_expr

) EXIT ...

END DO !

!

! Code block The block of statements between the DO and END DO are repeated indefinitely until the

logical_expr

becomes true and the EXIT statement is executed. After the EXIT statement is executed, control transfers to the first statement after the END DO.

Statistical Analysis Example: The average and standard deviation of set of numbers ( xi is sample i out of N samples):

Problem 1. State the problem Calculate the average and the

standard deviation of a set of , if measurements are either positive or zero, and if we do not know in advance how many measurements there are. A negative input value will mark the end of the set of measurements.

2. Define the inputs and outputs The inputs are a number of positive or zero real (floating-point) numbers.

Algorithm 3.

Design the algorithm

(steps) • Accumulate the input data • Calculate the mean and standard deviation • Write out the mean, standard deviation, and number of points

PseudoCode • The pseudocode for these steps is: • Initialize n, sum_x, and sum_x2 to 0 • WHILE • • Prompt user for next number Read in next x • • IF x < 0. EXIT • n  n + 1 sum_x  sum_x + x • sum_x2  sum_x2 + x**2 • End of WHILE

PseudoCode • Pseudocode to Calculate the mean and standard deviation. • x_bar  sum_x / REAL(n) • std_dev  SQRT((REAL(n)*sum_x2 - sum_x**2) / (REAL(n)*REAL(n-1))) • Write out the results : • Write out the mean value x_bar • Write out the standard deviation std_dev • Write out the number of input the n data points

Turn the algorithm into Fortran statements

PROGRAM stats_1 IMPLICIT NONE INTEGER :: n = 0 ! The number of input samples.

REAL :: std_dev = 0. ! The standard deviation of the input samples.

REAL :: sum_x = 0. ! The sum of the input values. REAL :: sum_x2 = 0. ! The sum of the squares of the input values. REAL :: x = 0. ! An input data value.

REAL :: x_bar ! The average of the input samples .

Turn the algorithm into Fortran statements

DO ! While Loop to read input values.

! Read in next value WRITE (*,*) 'Enter number: ' READ (*,*) x WRITE (*,*) 'The number is ', x ! Test for loop exit IF ( x < 0 ) EXIT ! Otherwise, accumulate sums.

n = n + 1 sum_x = sum_x + x sum_x2 = sum_x2 + x**2 END DO

Turn the algorithm into Fortran statements

! Calculate the mean and standard deviation x_bar = sum_x / real(n) std_dev = sqrt( (real(n) * sum_x2 - sum_x**2) / (real(n) * real(n-1)) ) ! Tell user.

WRITE (*,*) 'The mean of this data set is:', x_bar WRITE (*,*) 'The standard deviation is: ', std_dev WRITE (*,*) 'The number of data points is:', n END PROGRAM stats_1

Test the program

For input values: 3, 4, and 5 :

3.

The number is 3.000000

Enter number:

4.

The number is 4.000000

Enter number:

5.

The number is 5.000000

Enter number:

-1.

The number is -1.000000

The mean of this data set is: 4.000000

Enter number: The standard deviation is: 1.000000

The number of data points is: 3

The DO WHILE Loop

There is an alternate form of the while loop in Fortran 95/2003, called the DO WHILE loop. The DO WHILE construct has the form DO WHILE (

logical_expr

) ...

! Statement 1 ...

...

...

END DO ! Statement 2 ! ...

! Statement n Good Programming Practice : Do not use DO WHILE loops in new programs. Use the more general while loop instead.

The Iterative or Counting Loop DO index = istart, iend, incr Statement 1 !

...

Statement n !

! Body END DO

Loop index variable

: index

Loop parameters (constants or variables) :

istart, iend, incr

Note: All these variables or constants must be declared before the loop as INTEGER .

The Iterative or Counting Loop How loop works: 1. Evaluate parameters: istart, iend, and incr 2. Index  istart; If index*incr <= iend*incr, the program executes the statements within the body of the loop. 3. recalculate: index = index + incr 4. If index*incr <= iend*incr, goto to 2. above and repeat .

The Iterative or Counting Loop

The number of times

the count loop (without an IF () EXIT)

is executed

is given by the integer arithmetic expression:

(iend-istart+incr)/incr

Also, the count loop

will not be executed

if

for all index values for which: index*incr > iend*incr

DO i = 1, 10 Statement 1 ...

Statement n END DO DO i = 1, 10, 2 Statement 1 ...

Statement n END DO Examples

DO i = 1, 10, -1 Statement 1 ...

Statement n END DO DO i = 3, -3, -2 Statement 1 ...

Statement n END DO Examples

Examples

The Factorial Function:

function is defined as: To illustrate the operation of a counting loop, we will use a DO loop to calculate the factorial function. For integer numbers n (in mathematics) the factorial n! =1 if n = 0 n! = n * (n-1) * (n-2) * ... * 3 * 2 * 1 if n >= 1 e.g. for n=4 then n!=4*3*2*1=24.

The Fortran code to calculate n factorial for n>0 is: n_factorial = 1 ! We use the variable name n_factorial for n!

DO i = 1, n n_factorial = n_factorial * i END DO

Examples If n is 5, the DO loop parameters will be istart = 1, iend = 5, and incr = 1. The factorial will be: 1*2*3*4*5= 120

Example 4-3—Calculating the Day of Year:

1. Years evenly divisible by 400 are leap years. 2. Years evenly divisible by 100 but

not

by 400 are not leap years. 3. All years divisible by 4 but

not

by 100 are leap years. 4. All other years are not leap years. See the program Figure 4-6 (

PROGRAM doy

).

Examples (Programming Pitfalls) DO i = 3, 2 ...

END DO Since: (iend-istart+incr)/incr= (2-3+1)/1=0, the DO loop will never be executed. Count down loops: DO i = 3, 1, -3 ...

END DO Since:(iend-istart+incr)/incr= (1-3+(-3) )/(-3) = -5/(-3) = 1, the DO loop will be executed only once (for i=3).

Examples (Program Pitfalls)

Good Programming Practice:

Never modify the value of a DO loop index variable while inside the loop. Never depend on a loop index variable to retain a specific value after a DO loop completes normally. So, avoid using it in statements after the loop (see the examples in the next page).

Examples (Program Pitfalls) INTEGER :: i ! Branch out of the loop before completion DO i = 1, 5 ...

IF (i >= 3) EXIT ...

END DO

The value of the index variable i=3 when the loop is completed.

INTEGER :: i DO i = 1, 5 ...

END DO WRITE (*,*) i

The value of the index variable i is undefined when the loop is completed.

The CYCLE and EXIT Statements • • There are two additional statements which can be used to control the operation of while loops and counting DO loops: CYCLE and EXIT. If the CYCLE statement is executed in the body of a DO loop, the execution of the

current iteration

of the loop will stop, and control will be returned to the top of the loop. The loop index will be incremented, and execution will resume again if the index has not reached its limit.

The CYCLE and EXIT Statements Example with CYCLE: PROGRAM test_cycle INTEGER :: i DO i = 1, 5 IF ( i == 3 ) CYCLE ! If true go next to start of the loop WRITE (*,*) i END DO WRITE (*,*) ‘End of loop!’ END PROGRAM test_cycle

The CYCLE and EXIT Statements If the EXIT statement is executed in the body of a loop, the execution of the loop will stop and control will be transferred to the first executable statement after the loop. An example of the EXIT statement in a DO loop is shown below . PROGRAM test_exit INTEGER :: i DO i = 1, 5 IF ( i == 3 ) EXIT ! If true go next to end of the loop WRITE (*,*) i END DO WRITE (*,*) ‘End of loop!’ END PROGRAM test_exit

Named Loops

[name:] DO Statement Statement ...

IF ( logical_expr ) CYCLE [name] ...

IF ( logical_expr ) EXIT [name] ...

END DO [name]

Good Programming Practice:

Assign a name to any large and complicated loops in your program to help you keep the parts of the construct associated together in your own mind

.

Named Loops

[name:]

DO index = istart, iend, incr Statement Statement ...

IF (

logical_expr

) CYCLE

[name]

...

END DO

[name]

Nesting Loops It is possible for one loop to be completely inside another loop. If one loop is completely inside another one, the two loops are called

nested loops

. The following example shows two nested DO loops used to calculate and write out the product of two integers.

Nesting Loops PROGRAM nested_loops INTEGER :: i, j, product DO i = 1, 3 DO j = 1, 3 product = i * j WRITE (*,*) i, ' * ', j, ' = ', product END DO END DO END PROGRAM nested_loops

Nesting Loops The results are 1 * 1 = 1 1 * 2 = 2 1 * 3 = 3 2 * 1 = 2 2 * 2 = 4 2 * 3 = 6 3 * 1 = 3 3 * 2 = 6 3 * 3 = 9

When a Fortran compiler encounters an statement with the innermost currently open loop. Good Programming Practice:

understand and debug. END DO

statement, it associates that

Assign names to all nested loops so that they will be easier to

Nesting Loops Compiling a named loop helps us find errors PROGRAM bad_nested_loops_2 INTEGER :: i, j, product outer: DO i = 1, 3 inner: DO j = 1, 3 product = i * j WRITE (*,*) i, ' * ', j, ' = ', product END DO outer END PROGRAM bad_nested_loops_2 bad_nested_loops_2.f90(7) : Error: The block construct names must match, and they do not. [OUTER] END DO outer -------^ bad_nested_loops_2.f90(3) : Error: An unterminated block exists.

outer: DO i = 1, 3 ^ compilation aborted for bad_nested_loops_2.f90 (code 1)

Nesting Loops If two loops are to be nested, one of them must lie completely within the other one. e.g. of DO loops are incorrectly nested outer: DO i = 1, 3 ... inner: DO j = 1, 3 ... END DO outer ... END DO inner Good Programming Practice: Assign names to all nested loops so that they will be easier to understand and debug.

Nested Loops with CYCLE or EXIT If a CYCLE or EXIT statement appears inside an nested loops, then the CYCLE or EXIT statement refers to the

innermost

of the loops containing it. e.g.

unnamed

set of PROGRAM test_cycle_1 INTEGER :: i, j, product DO i = 1, 3 DO j = 1, 3 IF ( j == 2) CYCLE product = i * j WRITE (*,*) i, ' * ', j, ' = ', product END DO END DO END PROGRAM test_cycle_1

Nested Loops with CYCLE or EXIT The resulting output values are 1 * 1 = 1 1 * 3 = 3 2 * 1 = 2 2 * 3 = 6 3 * 1 = 3 3 * 3 = 9

Nested Loops with CYCLE or EXIT PROGRAM test_cycle_2 ! Skip remainder of the code block of the

inner

DO INTEGER :: i, j, product outer: DO i = 1, 3 inner: DO j = 1, 3 IF ( j == 2) CYCLE outer product = i * j WRITE (*,*) i, ' * ', j, ' = ', product END DO inner END DO outer END PROGRAM test_cycle_2 The resulting output values are 1 * 1 = 1 2 * 1 = 2 3 * 1 = 3

Nesting Loops within IF Constructs • • Nesting Loops within IF Constructs and vice versa It is possible to nest loops within block IF constructs or block IF constructs within loops. If a loop is nested within a block IF construct,

the loop must lie entirely within a single code block

of the IF construct.

Nesting Loops within IF Constructs e.g. illegal or incorrect nesting: outer: IF ( a < b ) THEN ... inner: DO i = 1, 3 ... ELSE ... END DO inner ... END IF outer

Nesting Loops within IF Constructs e.g. Correct nesting: outer: IF ( a < b ) THEN ... inner: DO i = 1, 3 ... END DO inner ... ELSE ... END IF outer

Character assignments and substrings e.g. CHARACTER(len=3) :: file_ext file_ext = 'f' ! Stores ‘f ‘ into variable file_ext e.g. CHARACTER(len=3) :: file_ext_2 file_extent_2 = 'FILE01' ! Stores ‘FIL‘ into var. file_ext e.g. str1 = '123456' then substring str1(2:4) is '234'

Example 4-5

What will the contents of variables a, b, and c be at the end of the following program?

Character assignments and substrings PROGRAM test_char1 CHARACTER(len=8) :: a, b, c a = 'ABCDEFGHIJ' ! a stores 'ABCDEFGH' b = '12345678' ! '12345678' c = a(5:7) ! c stores 'EFG ' (it has 5 blanks) b(7:8) = a(2:6) ! b stores '123456BC' WRITE(*,*) " a =", a WRITE(*,*) " b =", b WRITE(*,*) " c =", c END PROGRAM test_char1

Concatenation operator (//) PROGRAM test_char2 CHARACTER(len=10) :: a CHARACTER(len=8) :: b, c a= 'ABCDEFGHIJ' b = '12345678' c = a(1:3) // b(4:5) // a(6:8) ! c stores 'ABC45FGH' END PROGRAM test_char2

Relational Operators with Character Data Character strings can be compared in logical expressions using: ==, /=, <, <=, >, and >= : e.g. '123' == '123' is true e.g. '123' == '1234' is false Order is the

collating sequence

of the characters on the computer: It is the order number in the ASCII set of symbols: e.g. 'A' is 65 and 'B' is 66; also 'a' is 97, and ‘b' is 98 Comparison is left to right and character per character: e.g. true: 'AAAAAB' > 'AAAAAA' 'AB' > 'AAAA' and 'AAAAA' > 'AAAA'

Character intrinsic functions Book Table 4-1 character functions: ACHAR(ival) IACHAR(char), LEN(str1), LEN_TRIM(str1) , TRIM(str1) See also appendix (B.7 pp. 899-902).

Skip program Example 4-6—Shifting Strings to Upper Case

. One can compare character strings using the intrinsic functions: LGE(str1,str2), LGT(str1,str2), LLE(str1,str2),LLT(str1,str2) e.g. Fortran statements: str1='AAAAAB' , str2= 'AAAAAA' WRITE(*,*) LGT(str1,str2) ! It prints TRUE

Character intrinsic functions e.g.Quiz 4-2 What is written out by each of the WRITE statements below? PROGRAM test_char CHARACTER(len=10) :: str1 = 'Hello' CHARACTER(len=10) :: str2 = 'World' CHARACTER(len=20) :: str3 str3 = str1 // str2 WRITE (*,*) LEN(str3) ! 20 WRITE (*,*) LEN_TRIM(str3) ! 15 str3 = TRIM(str1) // TRIM(str2) WRITE (*,*) LEN(str3) ! 20 WRITE (*,*) LEN_TRIM(str3) ! 10 END PROGRAM test_char

Example 4-7—Physics--The Flight of a Ball

Range of ball thrown in the air will be: range= -2.0*v0*v0 /g*cos(theta)*sin(theta) Write a program with a loop to check angles theta = from 1 to 90 degrees PseudoCode: DO for theta = 0 to 90 degrees Convert theta to radians Calculate the range of the ball for each angle theta Determine if this theta yields the maximum range so far Write out the range as a function of theta END of DO WRITE out the theta yielding maximum range See (FIGURE 4-13) PROGRAM ball

Chapter 5

Basic I/O Concepts

Objectives

• Formatted WRITE • Format descriptors : I, F, E, ES, L, A, X, T, / • Formatted READ • Files

5.3 Format descriptors

Table 5-2 : Symbols with format descriptors Symbol meaning c column number d no. of digits to the right of the decimal place of the real input/output m minimum no. of digits to be displayed n no. of spaces to skip r Repeat count (i.e. no. of times to use a group of descriptors) w Field width (i.e. no. of characters to use in input or output)

Integer output : rIw e.g. INTEGER :: t,d,h

t= 23 d= 715 h= -12 WRITE(*,200) t,d,h 200 FORMAT(I5,I5,I5) will produce ♦♦♦ 23 ♦♦ 715 ♦♦ -12

Real output : rFw.d

e.g.

x = 3.14159

y = -275.3024

z = 12.9999

WRITE(*,200) x, y, z 200 FORMAT(F10.3, F10.3, F10.3) !using r=3: 200 FORMAT(3F10.3)

will produce ♦♦♦♦♦

3.142

♦♦

-275.302

♦♦♦♦

13.000

Real output : rEw.d

The

E

edit descriptor, produces a representation of a real number consisting of a decimal mantissa,

fr

, in the range 0.1 <

fr

< 1.0, with

d

digits of

fr

, followed by a four character exponent e.g.

For d=4 : + 0.d1d2d3d4E+ee

the whole number will occupy a field width of

w

characters. So, w >=d+7 (

7 characters

:

mantissa sign, 0, . , E, exp.sign, 2 for exp.

)

Real output : Ew.d

e.g. The number -0.000036176 4, will be output as:

E10.4

E12.6

-

.3618E-04 .361764E-04 E14.8

-

.36176400E-04

But with Fw.d format, we lose digits in the display:

F10.4

♦♦♦ -

0.0000

F12.6

F14.8

♦♦♦ -

0.000036

♦♦♦ -

0.00003618

Real output : ESw.d

The

ES

edit descriptor, produces a representation of a real number consisting of a decimal mantissa,

fr

, in the range 1.0<

fr

< 10.0, with

d

digits of

fr

, followed by a four character exponent e.g.

For d=3: + d1.d2d3d4E+ee

the whole number will occupy a field width of

w

characters. Again, w >=1+d+6 (6 characters:mantissa sign, . , E, exp.sign, 2 for exp.)

Real output : ESw.d

e.g. The number -0.0000361764, will be output as:

ES10.3

ES12.5

ES14.7

-

3.618E-05 3.61764E-05 3.6176400E-05

Logical output: rLw

Output

w

- 1 blanks, followed by

T

or

F

to represent a logical value e.g.

LOGICAL :: output =.TRUE. ,debug=.FALSE.

WRITE(*,’2L5’) output, debug

will produce: ♦♦♦♦

T

♦♦♦♦

F

Character output: rA or rAw

Aw

A

Output a character string in the next w character positions.

Output a character string, starting at the next character position, with no leading or trailing blanks.

e.g:

CHARACTER(LEN=26) :: long_name = "Llanfair...gogogoch" WRITE(*,”(A19)”) long_name produces:

Llanfair...gogogoch

Horizontal positioning: X or T

nX Ignore the next

n

character positions.

Tc position

c

Output the next item starting at character e.g.

REAL :: x=1.,y=2.,z=3.

WRITE(*,200) "X = ", x, " Y = ", y, " Z = ", z 200 FORMAT (A4,F4.1,2(2X,A5,F4.1)) will produce:

X= 1.0 Y= 2.0 Z = 3.0

Example of making a table using a do loop

An example of tabular printing:

PROGRAM tabular_output IMPLICIT NONE REAL, PARAMETER :: third=1.0/3.0

REAL :: x INTEGER :: i WRITE(* ‘The output….’ DO i=1,10 ! WRITE(* ‘The output….’ x=i WRITE(*, '(F15.4,3X,F15.4,3X,F15.4)‘) x, SQRT(x), x**third END DO END PROGRAM tabular-output EXAMPLE: Fig 5-5 (PROGRAM table)

Example of making a table using a do loop

The output from this program will be

1.0000

2.0000

1.0000

1.4142

1.0000

1.2599

3.0000

4.0000

5.0000

6.0000

7.0000

1.7321

2.0000

2.2361

2.4495

2.6458

1.4422

1.5874

1.7100

1.8171

1.9129

8.0000

9.0000

10.0000

2.8284

3.0000

3.1623

2.0000

2.0801

2.1544

Example of ES format EXAMPLES: Fig 5-8 (PROGRAM capacitor)

Equations type: Either (1) or (2) : (1) Given capacitance in farads and voltage in volts compute (total) charge in coulombs: Q = C*V (2) C= Q/V Also, compute energy in joules: E = (1/2)*C*V*V And no. of electrons (6.241461*10^(18) ) per coulomb of charge.

e.g. Test with input data: Volt.=100.00 V, Total charge=0.01 C (i.e. type =2 case), Output is: Capacitance = 1.000E-04 F, No. of electrons=6.241E+16 Total energy =.5000 joules.

READ Format

Edit descriptors

for input:

Descriptor Meaning

I

w

F

w

.d

[ Read the next

w

characters as an integer (e.g. w=8) Read the next

w

characters as a real number with

E A L

n w.

Aw w

X

Tc

d d

digits after the decimal place (e.g. w=14 , d=5) if no decimal point is present] Read the next

w

characters as characters Read sufficient characters to fill the input list item, stored as characters Read the next

w

characters as the representation of a logical value Ignore the next

n

characters (e.g. n=2) Next character to be read is at position

c (e.g. c=5)

Examples of READ Formats

Format I

w: Assume that we type on the keyboard 123456789 (as input):

READ (*, '(4X, I5)‘) num !

Leave 4 spaces then read num=56789 will ignore the first four columns and then read;

READ (*,'(I2, 3X, I3)‘) i, j !

Read i, leave 3 spaces; read j, will cause the value 12 to be stored in i and 678 in j.

Examples of READ Formats Assume that we type on the keyboard 123456789 (as input):

READ (*,'(T4, I2, T8, I2, T2, I4)‘) x, y, z !

Move to pos. 4, read x,… will read the number 45 into

x

, reading 89 into

y

, and reading the number 2345 into

z

.

Examples of READ Formats

Format F

w.d: The effect of the

F

edit descriptor during input with decimal point present.

READ (*,'(3F10.4)') a,b,c

1.5 0.15E+01 15.0E-01 ---------|---------|----------|---------| 10 20 30 40 Then a , b, c store the value 1.5

Examples of READ Formats The effect of the

Fw.d

edit descriptor during input without decimal point present. Then the d last digits are the fraction part.

READ (*,'(3F10.4)') a,b,c

15 150 15000 ---------|---------|---------|---------| 10 20 30 40 Then : a is 0.0015, b is 0.0150, c is 1.5000

Explanation

: Due to w=10, the first 10 places will be read for ‘a’ and due to d=4, there will be 4 fractional digits , similarly for b ,c

Examples of READ Formats

Good Programming Practice

: Always include a decimal point with REAL inputs, when using formated READ.

Format L

w: The edit descriptor used with logical data, takes the form: Lw (e.g. w=4) This edit descriptor processes the next

w

characters to derive either a

true

value, a

false

value, if the 1 st character is T or F. Else an error will occur.

Examples of READ Formats

Format A

: An

A

edit descriptor without any field width

w

treated as is though the field width was identical to the length of the corresponding input list item. Thus, for

CHARACTER (len=10) :: ch1 CHARACTER (len=8) :: ch2 CHARACTER (len=15) :: ch3

the following two statements will have an identical effect:

READ(*, '(A10, A8, A15)‘) ch1, ch2, ch3 READ (*,'(A, A, A)',) ch1, ch2, ch3 eg. Inputline:Mississipi Amazon longest river

5.5 Files and File processing

The UNIT

The UNIT number indicates which "stream" of input or output you desire to access.

By default the standard streams are keyboard (input) and monitor (output).

You can change the output to stream to (or input to stream from) a file of choice.

The number of the UNIT appears to be up to you to decide.

The FILE

This is the name of the file you wish to open.

You may also reference a character variable that holds the name of the file to access.

We call ‘Record’ one line of data the file.

5.5 Files and File processing

The STATUS

The STATUS specifies the status of the file to be opened.

The status of the file may be: – OLD - file exists – NEW - file needs to be created – REPLACE - file needs to be essentially deleted and recreated – SCRATCH - a temporary file that will be deleted when execution completes

The ACTION

Why you open the file? READ, WRITE, or READWRITE?

The IOSTAT tells

if all was correct or an error occur

5.5 Files and File processing

The OPEN Statement

OPEN(UNIT = 8, FILE = "data", STATUS = "OLD", ACTION = "READ", IOSTAT = ierror) The statement above will open the file named

data

reading.

for The variable

ierror

is an integer that will hold the status of the open the request. If open was successful.

ierror

equals zero then the UNIT is set to 8. this number is arbitrary, but it is unique among open I/O streams.

5.5 Files and File processing Examples of WRITE and READ with files: Figures 5-10 (PROGRAM

read_file

), 5-11 (PROGRAM

scratch_file

), 5-12 (PROGRAM

least_squares_fit

)

Figure 5.12 Least-squares fit

DATA FITING BY LEAST SQUARES :

Assume that we have a (set) table of data from experiments e.g and we want to pass a line as near to the data as possible. If we know that the data satisfy : y = ax + b; (

Note in our book: ‘a’ is ‘m’ and ‘b’ is the same as ‘b’ here)

TABLE x y x 1 x 2 y y 1 . .

2 x n y n

We know that if we have n = 2 pairs of data i.e.

then there is one line connecting them (if x 1 != x 2 ) TABLE x y x 1 x 2 y y 1 2 BUT if we have n > 2, then it is not possible (expect special cases) for the line to connect the (x,y) points. Then we want a line to pass ‘as near as possible’ y y7 y1 0 x1 x2……..x7 x

The problem of finding a line that passes as near as possible is solved by the ‘least squares’ method which finds the a,b which make MINIMUM the sum [ax i + b – y i ] 2 From LINEAR ALGEBRA and CALCULUS we know that the Solution (assuming two x i ’s at least are different) is

a b

  [(( 

n

((

i

 

n i

1  1

y i x i

)  )(

i a



n y

 ( 1

i



n

 1

i

) 

n

(

x i

)) /

n i



n

 1

x y i i

)) /(( 

n

 1

i x i

) 2 

n

(

i



n

 1

x i

2 ))]

Note: solve for ‘a’ first and use it in ‘b’. Note in our book: ‘a’ is ‘m’ and ‘b’ is the same as ‘b’ here.

Program in Figure 5-12.

For each (x i , y i ) we call ‘residual’ r(x i ) = [ax i + b – y i ] ; this states how ‘closely’ the line y = ax + b passes through the point (x i , y i ) e.g. y 10*r(x1) 0 x x1 x2 10*r(x2) The ‘residual sum’:

i n

  1 x7 (

i

) 2 

i n

  1 [

ax i



b



y i

] 2 gives the ‘goodness of the fit’. For best fit the residual sum

i



n

 1

r x i

2  0

Chapter 6 Arrays

Motivation:

We know REAL or INTEGER variables or constants

e.g.

REAL :: x REAL, PARAMETER ::

pi

=3.14

• These allow us to program computations with reals or integers.

• How about working with vectors or matrices (needed in linear algebra)?

Chapter 6 Arrays

e.g: X 3x1 or = [X 1 A 3x1 = [a 11 X 2 a 12 X 3 ] T = [1.0 2.0 3.0] T a 13 ] T A 2x2 = [a 11 a 12 a 21 a 22 [1.0 0.0 ] = 2.0 4.0] = a(1:2, 1:2) A 2x3 = [a 11 a 21 a 12 a 22 a 13 [0.5 0 0.3 a 23 ] = 0.1 1.0 1.2] = a(1:2, 1:3)

The array concept

So far we have used one name to refer to one location in the computer's memory.

A group, or array

, of locations in the memory, are identified by the same name but with an

index, or subscript

, to identify individual locations.

As the vector A, the individual elements are A1, A2, ... An We call such an ordered set of related variables an array, and we refer to the individual items within the array as array elements.

The array concept

We follow the name of the array by an identifying integer value enclosed in parentheses: A(1), A(2), ..., A(n) An array element is defined by writing the name of the array followed by a subscript, where the subscript consists of an integer.

Expression (known as the subscript expression) is enclosed in parentheses. Thus, x(10) y (2*i+4)

The array concept

Array (in programs) is a name of a variable or constant with an index set . e.g. [a 1 , a 2 , a 3 ] (of linear algebra) is written as : a(1),a(2), a(3) (in FORTRAN 90).

More generally a(1), a(2),….a(n) where n is an integer value (e.g. n=3 for previous example) An array element is (a subscript expression) e.g: x(3), y(2* i +4)

Note:

subscript expression must be integer ; but it can be an array itself or contain an array as part of it.

Chapter 6 Arrays

Declaring Arrays

Before an array can be used, its

type

and the

number of elements

it contains must be declared to the compiler.This must be done so that the compiler knows what type of data are to be stored in the array and how much memory is required to hold the indicated number of elements of that type.

Example from the book

can hold 16 elements.

. Let's assume we need to store the voltage (v) readings from 16 different experiments, we could declare 16 individual REAL variables or one REAL array that REAL :: v1, v2, v3, v4, v5, v6, v7, v8, v9, v10 REAL :: v11, v12, v13, v14, v15, v16 Probably not a good idea. Try the array instead! REAL, DIMENSION(16) :: voltage

Chapter 6 Arrays • Each individual REAL value is represented by an index (or subscript).

Index: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

Value: 3.2 2.4 1.0 1.1 3.2 2.2 1.2 4.3 2.3 4.5 3.2 1.0 0.3 3.4 3.4 8.3

6.1 Array declarations The declaration of scalar (e.g. REAL) variables causes the compiler to allocate an appropriate storage unit to contain its value. When we declare an array variable, the compiler will need to allocate several storage units.

There are two ways using either a

dimension attribute

or an

array specification

applied to the variable name.

Consider e.g.

REAL, DIMENSION(50) :: a, b, c REAL :: a(50), b(50), c(50) REAL, DIMENSION(50) :: a,b,c, x(20),y(20),z

To store arrays in memory we need: - name (

e.g.

a) - type (

e.g.

REAL) - Dimension i.e number of elements e.g. ‘3’ for a(1), a(2), a(3).

FORTRAN declaration for arrays: e.g. REAL,DIMENSION(50) :: a,b,c 50 means: (1:50) -Mixed dimensions: REAL,DIMENSION(50) :: a,b,c,x(20),y(20),z or better

REAL, DIMENSION(50) :: a,b,c,z REAL, DIMENSION(20) :: x,y

6.1

NOTE:

This means e.g. we have stored x(1),x(2),x(3),………x(20) i.e. lb :lower bound=1 y(1),y(2),y(3),………y(20) ub : upper bound=20 By default the subscripts will start at 1, but if we wish the subscripts to have a different range of values then we may provide the

lower bound (lb)

and the

upper bound (ub)

explicitly, separated by a colon:

(lb:ub) REAL, DIMENSION(11:60) :: a,b,c ![extent=no of elements REAL, DIMENSION(-20:-1) :: x ! = ub-lb+1] REAL, DIMENSION(-9:10) :: y REAL, DIMENSION(0:49) :: z

6.1 Five technical terms that are of great importance when discussing arrays in Fortran.

 Fortran permits up to seven subscripts, each of which relates to one

dimension

of the array. For each dimension, there are two bounds, the lower bound and the upper bound.

 The number of permissible subscripts for a particular array is called its

rank.



The extent of a dimension is the number of elements in that dimension.



The size of the array is the total number of elements

.

 The

shape

of an array is determined by its rank and the extent of each dimension.

I

Chapter 6 Arrays

6.2.1 Array Elements Are Just Ordinary Variables

Each element of an array is a variable just like any other variable! e.g. program segment …… INTEGER, DIMENSION(10) :: arrayInt REAL, DIMENSION(3) :: arrayReal arrayInt(1) = 5 arrayReal (3) = arrayInt(1)/4.0

WRITE(*,*) "arrayInt(1) = ", arrayInt(1) WRITE(*,*) "arrayReal(3) = ", arrayReal(3)

Initialization of Array elements Four different approaches for initialization: (1)

Using a loop

: (e.g. program segment) INTEGER, DIMENSION(10) :: arrayInt INTEGER :: lcv DO lcv = 1, 10, 1 arrayInt(lcv) = lcv * 2 END DO

e.g. PROGRAM squares ! Fig 6-2 in book

Initialization of Array elements (2)

Provide values for each of these elements by means of an array constructor.

An array constructor consists of a list of values enclosed between special delimiters, (/ …. /):

(/ value_l, value_2, … /)

If

arr

is an integer array of size 10, its elements could be set to

arr = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 /)

Chapter 6 Arrays (3)

Use implied DO loop

: e.g.: Create an integer array (dimension 1000) and initialize each element to be equivalent to its index. e.g.

INTEGER, DIMENSION(1000) :: arrayInt = (/ (i, i=1, 1000) /)

e.g. PROGRAM square_roots ! Fig 6-3 in book Another Example (text page 256):

INTEGER, DIMENSION(25)::array4=(/ ( (0,i=1,4), 5*j, j=1,5 ) /) The inner loop (0,i=1,4) executes completely for each outer loop Index j. So the resulting array4 is: 0,0,0,0,5, 0,0,0,0,10, 0,0,0,0,15, 0,0,0,0,20, 0,0,0,0,25

Chapter 6 Arrays (4)

Use whole array operation

: Initialize an array (dimension 5000) to zero for all elements.

INTEGER, DIMENSION(5000) :: arrayInt =1.0

! OR: INTEGER, DIMENSION(5000) :: arrayInt arrayInt = 1.0

Chapter 6

(

Changing the Subscript Range of an Array)

We can change the range of values that an array's subscript can take when we declare the array.

INTEGER, DIMENSION(-50:49) :: arrayInt !! Or INTEGER:: lower_bound = -50, upper_bound=49 INTEGER, DIMENSION(lower_bound:upper_bound) :: arrayInt

Shape

of an array is the number of dimensions (here it is 1) and the extent of each dimension. We have only seen one dimensional arrays so far.

Extent

of an array is the upper_bound - lower_bound + 1

e.g. PROGRAM squares_2 ! Fig 6-4 in book

Chapter 6 Arrays (out of bounds subscripts) The Fortran compiler may not report ‘bounds checking’errors at compile time. Then we may get ‘unwanted’ results.

Example

:

PROGRAM bounds! Fig 6-6 in book

REAL, DIMENSION(5) :: a This declaration means that the bounds are: 1:5 for array Thus we can not read/write other elements e.g. a(6) does not really exist.

Use of Named Constants with Array Declarations

Instead of using a hard-coded value as the size of an array, it is more common to use a constant to represent the array size.

INTEGER, PARAMETER :: nstud= 100 INTEGER, DIMENSION(nstud) :: stuID INTEGER, DIMENSION(nstud) :: course

Example

:

PROGRAM extremes! Fig 6-7 in book

Whole Array Operations

If the two arrays are the same shape, then they can be used in ordinary arithmetic operations and the operation will be applied on an element-by-element basis (

whole arrays

).

Example

:

PROGRAM add_arrays ! Fig 6-9 in book

INTEGER :: i REAL :: a(4) = (/1, 2, 3, 4/) REAL :: b(4) = (/5, 6, 7, 8/) REAL :: c(4) REAL :: d(4)

Whole Array Operations

!!ADDING ELEMENT BY ELEMENT DO i = 1, 4, 1 c(i) = a(i) + b(i) END DO !!WHOLE ARRAY ADDITION d = a + b DO i = 1, 4, 1 WRITE(*,*) "c(",i,") = ", c(i), "d(",i,") = ", d(i) END DO Other examples:

Whole Array Operations

Other examples: One must first check that the arrays are of the same shape (

conformable

)

REAL, DIMENSION(20) :: a, b, c REAL :: d(10:29) a = 1.55

b = c*d ! Whole array processing !! or, the same computation using a loop DO i=1, 20 ! DO loop array operations b(i) = c(i)*d(i+9) ! difficult!

ENDDO

Whole Array Operations

Other examples e.g.

1)

whole array processing:

REAL :: a(1:20), b(0:19), c(10:29), d(-9:10) a = c*d ! Fortran 90 style !! Or, The same computation using a loop DO i=1,20 ! FORTRAN

77

style b(i-1) = c(i+9)*d(i-10) ! b= c*d END DO 2) REAL:: array_1(1:20), array_2(-9:10) array_1 = 10*array_2 3) REAL:: a=0.0, b=0.0

REAL, DIMENSION(50) :: c = 0.0, d = 0.0

Whole Array Operations Intrinsic functions (p.48):

Syntax with variables (i.e. not arrays): e.g.

b = SQRT(a) If an intrinsic function (p.48)

has an array as an argument then the result of the function reference will be an array with the same shape as the argument. But, we need conformable arrays:

e.g. 1) REAL:: array_1(1:20), array_2(-9:10) array_1 = SIN(array_2)

Thus, the statement

array_1 = SIN(array_2)

assigns the sine of each element of the array

array_2

to each corresponding element of the array

array_1

.

Note:

that the assignment starts from the lower index and goes to the upper and it is the same as a loop. e.g.

DO i = 1, 20 array_1(i) = SIN(array_2(i - 10)) END DO

Whole Array Operations

Also, an intrinsic function has more than one argument then they must all be conformable.

e.g. 2) REAL :: a(1:20), b(0:19), c(10:29), d(-9:10) REAL::e(1:20), arr_max(1:20) arr_max = MAX(100.0, a, b, c, d, e)

Thus the statement

arr_max = MAX(100.0, a, b, c, d, e)

will assign to the elements of

arr_max

the maximum value of the corresponding elements of the arrays

a

,

b, c, d

and

e,

or 100.0.

Whole Array Operations

Can we use comparisons with whole arrays? IF (a < b) THEN WRITE(*,*) "a < b" ELSE WRITE(*,*) "a >= b" END IF

Not possible ! Syntax Error!

Arrays Subsets (or Sections) These are specified by subscript-triplets for each dimension of the array. The general form is: [<

subscritpt1

>]:[<

subscritpt2

>][:<

stride

>] The section starts at <

subscritpt1

> and ends at or before <

subscritpt2

> <

stride

> is the

incremen

t (or

jump

) by which the elements are selected.

<

subscritpt1

>, <

subscritpt2

> and <

stride

> must all be scalar integer expressions.

e.g. INTEGER :: i=3, j=7 INTEGER, DIMENSION(10) :: a =(/1,2,3,4,5,-6,7,8,9,10/) a(:) is the entire array

a

a(i:j) is the array subset: [3,4,5,-6,7]

Arrays Subsets (or Sections) Other examples for: INTEGER :: i=3, j=7 INTEGER, DIMENSION(10) :: a =(/1,2,3,4,5, -6,7,8,9,10/) a(i:j:k) , for k=3, is the array subset [3,-6] a(8:3:-1) ,i.e. from a(8) to a(3) in steps of –1, is [3,4,5,-6,7,8] a(8:3) ,i.e. from a(8) to a(3) in steps of 1, is an empty subset a(i:) , from a(i) to upper_bound of

a,

i.e.: [3,4,5,-6,7,8,9,10] a(:j) , from lower_bound of

a

to a(j), i.e.: [1,2,3,4,5,-6,7] a(::3), from lower_bound of

a

to upper_bound of

a

with stride= 3 i.e.: [1,4,7,10] a(i:i), i.e. :[3]

Good Programming Practices with Arrays

• Are arrays necessary? Yes, e.g. for Linear Algebra vectors and matrices and Multi dimensional Calculus • Array dimensions must be declared using named constants. This leads to fewer errors!

• Arrays should be initialized before use. The results of using uninitialized arrays are unpredictable.

• Stay within the array bounds.

Input and output of whole arrays and sections • Before using an array in prog’s must be able to input data in arrays and get output e.g.:

PROGRAM inputoutputData INTEGER :: first= 1, last= 100 ! Must declare all arrays READ (*,*) first,last,(arr(i), i=first,last) !Input !output

WRITE (*, *) (arr(i), i=1,99,2) ,arr(3),arr(4) WRITE (*, *) arr END PROGRAM inputoutputData

Input and output of whole arrays and sections • Array elements are treated in just the same way as scalar variables.

• An array name may appear in an input or output list, in which case it refers to the whole array.

• Input data to arrays and output results from arrays (next):

Input and output of whole arrays and sections Part of an array may be identified by use of an implied

DO

, the item takes the form

(object_list, implied_do_control)

For example:

WRITE (*, *) (arr(i), i=1, 99, 2), arr(3), arr(4)

Input and output of whole arrays and sections It is permitted for one or more of the controlling values for an implied

DO

in an input statement to be themselves input by the same statement:

READ (*,*) first, last, (arr(i), i=first, last)

This form of input statement must be used with care.

NOTE

: (1)We must be careful w/ defining ‘

first

’ and ‘

last

’: It must:

first < last

, (2) If ‘

arr

’ dimensions are

arr(low:upper)

then also must:

low <= first

and

last <= upper

.

Input and output of whole arrays and sections Example of input data into array: (A safe way of reading data with an implied

DO

)

PROGRAM array_input IMPLICIT NONE INTEGER, PARAMETER :: lower=-50, upper=50 INTEGER :: first, last, i REAL, DIMENSION(lower:upper) :: arr . . .

READ (*,*) first, last

Input and output of whole arrays and sections

IF (first>=lower .AND. last<=upper) THEN READ *, (arr(i), i=first, last) ELSE WRITE(*, *) 'Invalid array subscript specification! ' .

END IF .

END PROGRAM array_input Book Example: program array_io Fig. 6-11

EXAMPLES (sorting data e.g. numbers)

Sorting (i.e. ordering) by

straight selection

: Index Initial order After first exchange After second exchange After third exchange 1 2 3 4 5 6 7 8 7 1 8 4 6 3 5 2 1 7 8 4 6 3 5 2 1 2 8 4 6 3 5 7 1 2 3 4 6 8 5 7 After fourth exchange After fifth exchange After sixth exchange After seventh exchange 1 2 3 4 6 8 5 7 1 2 3 4 5 8 6 7 1 2 3 4 5 6 8 7 1 2 3 4 5 6 7 8

Sorting Sorting Data: Algorithm “

selection sort

”.

The list is divided into two sublists, sorted and unsorted, which we will separate. Initially, the entire array constitutes the unsorted list.

We go through the unsorted list and find the smallest element. Then we swap that element with the element at the beginning of the unsorted (or end of sorted list) data .

Sorting

What does it mean to swap? Answer

INTEGER :: x = 3, y= 2 INTEGER :: temp temp = x x = y y = temp

Sorting Once we find the smallest element in the unsorted list and move it to the front position. We then move one to the right.

This adds one datum to the sorted list and takes one datum away from the unsorted list.

Each time we complete such a swap we say we have completed a

sort pass

.

With n elements in an array we need n-1 sort passes to properly sort the data

Sorting PROGRAM selectionsort IMPLICIT NONE INTEGER:: SIZE,temp,i,j, smallest INTEGER :: data(50) INTEGER :: num_it = 0 READ(*,*) SIZE DO i = 1, SIZE, 1 READ(*,*) data(i) END DO

Sorting DO i = 1, SIZE-1, 1 smallest = i DO j = i+1, SIZE, 1 num_it = num_it + 1 IF (data(j) < data(smallest)) THEN smallest = j END IF END DO temp = data(i) data(i) = data(smallest) data(smallest) = temp END DO

Sorting DO i = 1, SIZE, 1 WRITE(*,*) "data(", i, ") = ", data(i) END DO WRITE(*,*) "Number of iterations: ", num_it END PROGRAM selectionsort Also, Book program for

Selection Sort: Figure 6-15:

program

sort1

Chapter 7 Introduction to procedures • •

Motivation: Top-Down Design

This is the idea that with any problem you should try to break down the problem into smaller easier to program tasks.

The code for these smaller tasks is isolated in to what is called a procedure (subroutine or function

)

.

Chapter 7 Introduction to procedures The book defines as follows: –

Subroutine

- a procedure that can return multiple results through calling arguments.

–

Function

- a procedure that returns only a single value that is used in the evaluation of an expression.

Arguments

- The data that goes into (or comes out from) a procedure

Chapter 7 Introduction to procedures

7.1 Subroutines:

A subroutine is a procedure that receives its input values and returns its results through an

argument list .

SUBROUTINE some_name (some argument list) ! DECLARATION SECTION ! EXECUTION SECTION RETURN END SUBROUTINE some_name Note: "variables" in the argument are called dummy arguments

( memory is not allocated for them, they act as placeholders for the actual arguments when the subroutine is called ).

Chapter 7 Introduction to procedures (1) The

subroutine

itself behaves as normal program and it has both a

declaration

and

execution

section.

(2) When the subroutine is

called

program

shifts

control of the from the calling program to the subroutine. (3) When the subroutine ends the control returns in the calling program in the statement after the

CALL e.g.



CALL some_name (actual argument list)

Chapter 7 Introduction to procedures

Example (book):

Fig. 7-1 and 7-2:

calc_hypotenuse

Need

calling program

(called

test driver

) to run the subroutine.

SUBROUTINE calc_hypotenuse ( side_1, side_2, hypotenuse )

Note:

3 dummy variables (first 2 are input and 3 rd output)

Chapter 7 Introduction to procedures • The variables being used as arguments must be declared (

REAL

,

INTEGER

,

CHARACTER

,

LOGICAL

,

arrays

, ) • The intent of each argument must be declared as well.

–

INTENT IN

- The value is coming in to the procedure and cannot be changed.

–

INTENT OUT

- The variable is coming in to the procedure and the result assigned to the variable will be "returned" from the procedure.

– The other (local) variables (i.e. not arguments) of the subroutine must be declared too. (e.g. REAL :: temp ) – The keyword

RETURN

is optional

Chapter 7 Introduction to procedures

Example (for INTENT): calc_hypotenuse

Fig. 7-1 and 7-2: • • •

7.1.2 The INTENT Attribute : INTENT (IN)

- Dummy argument is used only to pass input data to the subroutine.

INTENT (OUT)

- Dummy argument is used only to return results to the calling program.

INTENT (INOUT)

- Dummy argument is used both to pass input data to the subroutine and to return results to the calling program e.g. More Examples:

bad_call

(Fig7-5),

array2

(Fig7-6),

intent_test (

next slide

,

not in the book programs

)

Chapter 7 Introduction to procedures PROGRAM intent_test ! Example of INTENT(IN) IMPLICIT NONE INTEGER :: x = 1 INTEGER :: y = 2 CALL testing(x, y) WRITE(*,*) "The value of x = ", x, " and the value of y = ", y END PROGRAM intent_test SUBROUTINE testing(a, b) IMPLICIT NONE INTEGER, INTENT(IN) :: a INTEGER, INTENT(OUT) :: b ! Error because next we change variable a , which is input argument.

! it will work if we change IN to INOUT a = 3 b = 4 END SUBROUTINE testing

Chapter 7 Introduction to procedures

7.2 Sharing Data Using Modules:

• We have seen how data can be exchanged with subroutines as arguments.

• Another way of sharing data among programs, subroutines, and functions is by use of a module (or modules).

• A

module

is a unit that contains the definitions and initial values of data we wish to share between program units (functions, subroutines).

Example:

shared_data, test_module

(Fig 7-8, Fig 7-9)

Chapter 7 Introduction to procedures A module exists in order to make some or all of the entities declared than one program unit.

within it accessible to more One very important use of modules relates to

global accessibility

of variables, constants and derived type definitions.

A module allows a defined set of variables and/or constants to be made available to any program units which access them by means of an appropriate

USE

statement.

Chapter 7 Introduction to procedures

MODULE global_data IMPLICIT NONE SAVE ! Constant declarations REAL, PARAMETER :: pi=3.1415926, piby2=pi/2.0

! Variable declarations REAL :: global_1, global_2, global_3, global_4, global_5 END MODULE global_data Note: SAVE

should always be included in any module which declares any variables.

Chapter 7 Introduction to procedures The

USE

statement comes after the initial statement

(PROGPAM, SUBROUTINE

or

FUNCTION)

but

before any other statements

:

SUBROUTINE component USE global_data IMPLICIT NONE !Any additional specification statements !Executable statements END SUBROUTINE component Note:

All procedures which use a module that includes a

SAVE

statement are sharing the same copy of any variables, constants, type definitions or other entities accessed from by

USE

association.

Chapter 7 Introduction to procedures

7.3 Module Procedures

• A

module

can also be used to share procedures. The called procedure is said to have an

implicit interface

if we use

CALL subroutine_name

(because the

calling program

does not know the

types, intent,

etc of the arguments at

compile-time

). We can change the

implicit interface

to

explicit interface

by using modules.

Example:

bad_call2

(Fig. 7-12), which allows the compiler to discover the error of missing ‘IN’ keyword in the subroutine: ‘my_subs’ in program ‘bad_call’

Chapter 7 Introduction to procedures

7.4 FORTRAN Functions

A FORTRAN

function

is a

procedure

whose result is a

single number, logical value, character string or array

.

This result can be used to form a FORTRAN expression.

The expression may be on the right side of an assignment statement.

There are two types of functions, intrinsic and user-defined.

Intrinsic functions are those functions built into a FORTRAN language, such as SIN(x) or LOG(x).

See Appendix B.

Chapter 7 Introduction to procedures User-Defined functions are functions defined by programmers (not really users) to meet a specific need not addressed by the standard intrinsic functions.

General form of user-defined functions

REAL FUNCTION name(argument list) !Declaration

statements !executable statements name = arithmetic expression RETURN END FUNCTION name

Chapter 7 Introduction to procedures

Note:

(1) A function is (invoked) called by including it in an expression ( e.g.

name = expression)

or

WRITE (*,*) name )

(2) The parentheses around the argument list are required even if the list is blank.

(

3) The type of the function must be declared both in the function procedure and the calling programs.

Examples in book: (1)

quadf

(2)

sinc

(Fig7-16), (Fig 7-13),

test_sinc test_quadf

(Fig-7-17).

(Fig 7-14),

Chapter 7 Introduction to procedures Example:

REAL FUNCTION cube_root (x) IMPLICIT NONE ! Function to calculate the cubic root of a positive real number ! Dummy argument declaration REAL, INTENT (IN) :: x ! Local variable declaration REAL :: log_x log_x = LOG(x) cube_root = EXP (log_x/3.0) RETURN END FUNCTION cube_root

Chapter 7 Introduction to procedures

PROGRAM function_demo IMPLICIT NONE ! Variable declarations; actual arg’s REAL :: a=1. ,b=2. ,c=27. ,d=1.

!REAL :: cube_root ! Not enough REAL, EXTERNAL :: cube_root ! Need this a = b*cube_root (c) +d WRITE(*,*) a END PROGRAM function_demo

Chapter 7 Introduction to procedures

NOT IN EXAMS: 7.5.1 Passing user-defined functions as arguments:

We can have arguments (in procedures) which are

functions

. For this (

and always with user-defined functions

), we must declare (in the calling programs) the function as

EXTERNAL

e.g.

REAL, EXTERNAL :: func Example

:

ave_value

(Fig. 7-18) ,

test_ ave_value Example

:

sub_as_arguments

(Fig. 7-20) ,

test_sub_as_arguments

(Fig 7-21) (Fig 7-19)

NOT COVERED:

We can also have arguments (in procedures) which are

subroutines.

For this we must also declare (in the calling programs) the subroutine as

EXTERNAL.

Chapter 8 Additional Features of Arrays( Skip: 8.4-8.5) •

Multidimensional Arrays

• The arrays we have seen so far are 1-dimension ( rank=1) arrays. We often visualize such arrays as being in a single row.

• However, it might be more useful to think of the 1 dimension as being in one column based on how arrays are stored in memory.

• A 2-dimension array can be visualized as a table made up of several rows and columns. If the columns are a single dimension array then a table (2 dimensional array) can be viewed as a series of 1 dimensional arrays (columns).

Chapter 8 Additional Features of Arrays Why think of a 1-dimension array as a column? Because in physical memory, arrays are stored by column, such that all of column 1's contents reside consecutively then column 2's, column 3's, and so on. Such arrangement in memory is known as

column major order

. Not all programming languages will be column major order, many (e.g. C) are

row major order

.

In Chapter 6, we discussed rank = 1 (1-dimension or 1-D) arrays.

e.g. REAL, DIMENSION (11:60) :: a To solve problems in linear algebra, we need many dimensions (i.e. rank >= 2).

13.1 Matrices and two-dimensional arrays:

A 3 x 4 matrix A 3x4 : A = 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 A 1,4 A 2,4 A 3,4 e.g A = 2 3 4 3 4 5 4 5 6 5 6 7

The matrix A can be defined in FORTRAN 90: REAL, DIMENSION (3,4) :: a

Note:

in DIMENSION (3,4) the no. of rows is 1 st , then the no. of columns.

Other types of arrays are possible: e.g. LOGICAL, DIMENSION (10,4) :: b e.g. b 10x4 = b(1,1) b(1,2) b(1,3) b(1,4) ………………………………….

.………………………………… b(10,1) b(10,2) b(10,3) b(10,4)

and the values b(i,j) = .TRUE. Or .FALSE. for i = 1,…,10; j = 1,…,4.

The array elements (being ‘scalars’ i.e. single elements) can be used in expressions exactly as the simple type (e.g. REAL, e.g INTEGER, LOGICAL) variables or constants.

e.g for REAL, DIMENSION (3,4) :: a a(3,4) = 2.0*a(3,4) + 1.0

DO i = 1,4

! Doubles a(3,4) and adds 1 to it ! Replace row 1 of a by a(1,i) = a(3,i)

END DO DO i = 1,3

a(i,2) = a(i,1)

END DO

row 3 of a. Row 3 is unaltered ! Replace column 2 of a by column 1 of a. Column 1 is unaltered

Basic concepts for dim => 1 arrays :

Recall:

rank

= No. of dim’s and taken from ‘DIMENSION()’: e.g. REAL, DIMENSION(8) :: a !rank = 1 INTEGER, DIMENSION(3,10) ::b !rank = 2 array size =No. of array elements = product of extents = 3*10 = 30 Note : e.g. for b: DIMENSION(3 ,10 ) same as DIMENSION(1:3, 1:10 ) and the ‘extents’ are 3,10 Instead of the actual dim’s we can use ‘named constants’ e.g. INTEGER, PARAMETER :: s1=3, s2=10 … INTEGER, DIMENSION(s1,s2) ::b

‘shape’

of array : rank , extents of dim’s e.g.

shape of array ‘b’ above : rank=2 and extents = (3, 10) -e.g. arrays with different index bounds and the

same shape

e.g. Arrays a1, b1 below have the same shape as arrays : a, b above : REAL, DIMENSION(11:18) :: a1 INTEGER, DIMENSION(5:7, -10:-1) :: b1 -Array element order (in FORTRAN) is by columns :e.g.

REAL , DIMENSION(4,3) :: arr WRITE( * ,*) arr arr(1,1) , arr(2,1), arr(3,1), arr(4,1), arr(1,2), arr(2,2) arr(3,2) , arr(4,2), arr(1,3), arr(2,3), arr(3,3), arr(4,3)

arr(1,1) arr(1,2) arr(1,3) arr(4,1) arr(4,2) arr(4,3)

13.3 Array constructors : e.g

REAL , DIMENSION(0:49) :: arr = (/-1, (0 , i = 1,48),1 /) means: arr(0) = -1 arr(1) = 0 .

.

arr(48) = 0 arr(49) = 1

Since array constr. are 1 Dim we use ‘RESHAPE’ to apply to Dim > 1 arrays e.g.

a =

RESHAPE( (/1.0 ,2.0, 3.0, 4.0, 5.0,6.0/), (/2 , 3/) ) We could use nested implied ‘DO’ e.g. for 11 21 12 22 by the following declaration statement

INTEGER :: i, j REAL , DIMENSION(2,2) :: a = & RESHAPE( (/ ((10*i+j , i = 1,2), j =1,2) / ) , (/ 2,2 /) ) Which is the same as the following: INTEGER :: i , j REAL, DIMENSION(2,2) :: a DO j = 1,2 DO i = 1,2 a( i,j) = 10*i+j ENDDO ENDDO

Input/output with arrays :

The ways for input/output with arrays: (1) List individual elements (2) List of individual elements under an implied/regular DOloop (3) Use an array name (e.g ‘x’) with no subscripts in a whole array operation e.g. REAL, DIMENSION(50 , 8) :: x (x has 50 rows, 8 col’s) WRITE(*, ‘ (8F8.2)’) x

will print out the data in the following way x(1,1) x(2,1) x(3,1) x(4,1) x(5,1) x(6,1) x(7,1) x(8,1) x(9,1) x(10,1) x(11,1) x(12,1) x(13,1) x(14,1) x(15,1) x(16,1) . . . . . . . .

. . . . . . . .

. . . . . . . .

x(49,1) x(50,1) x(1,2) x(2,2) x(3,2) x(4,2) x(5,2) x(6,2) x(7,2) x(8,2) x(9,2) x(10,2) x(11,2) x(12,2) x(13,2) x(14,2) . . . . . . . .

. . . . . . . .

which is not at all what was wanted! It does not print the array as a matrix!

On the other hand, the implied Doloop statement WRITE(*, ‘(8F8.2)’) ((x(i,j) , j = 1,8 ), i = 1,50) will cause the results to be printed in the correct arrangement (like a matrix) x(1,1) x(1,2) x(1,3) x(1,4) x(1,5) x(1,6) x(1,7) x(1,8) x(2,1) x(2,2) x(2,3) x(2,4) x(2,5) x(2,6) x(2,7) x(2,8) . . . . . . . .

. . . . . . . .

. . . . . . . .

WHOLE ARRAY OPERATIONS

- Need Conformable arrays : i.e. arrays of same shape - e.g. (i)

REAL , DIMENSION (10) :: p,q REAL , DIMENSION (10 : 19) :: r

- Operations with conform. arrays do not need loops e.g. (i)

DO i = 1, 10 p(i) = q(i) + r(i + 9) END DO

same as whole array operation:

p = q + r

e.g. (ii)

REAL, DIMENSION(10, 21) :: x REAL, DIMENSION(0:9, -10:10) :: y REAL, DIMENSION(11:20, -20:0) :: z shape of x : rank = 2, extents (10 , 21) Also, shape of y and z : rank = 2, and extents = (10 , 21) Then the Fortran 90 whole array statement: x = y + z +10.0

has exactly the same effect as the following nested DO loop:

DO i = 1, 10 DO k = 1, 21 x(i,k) = y(i-1, k-11) + z(i+10, k-21) + 10.0

END DO END DO

Rules (from chapter 6) for working with whole array operations: (1) Two arrays are conformable, if they have the same shape (2) A scalar, including a constant, is conformable with any array (3) All intrinsic operations are defined between conformable arrays

e.g. (i) REAL :: x(4) = (/ 0.0, 1.1, 2.2, 3.3 /) REAL :: y(4) y = SIN(x) (ii) x= y+z +10.0 +EXP(y)+COS(z) The result is the same as if the function were applied to each element of the array. The shape of the resulting array must be the same as the input argument. All intrinsic functions e.g. ABS, SIN, COS, TAN, EXP, LOG, LOG10, MOD, and SQRT may also accept arrays as input.

Chapter 8 Additional Features of Arrays •

Array subsets

e.g. REAL, DIMENSION (1:5,1:5) :: a, b, c (i) a(:,1) is the 1 st column (ii) a(1,:) is the 1 st row (iii) a(1:3,1:5:2) takes an array subset of stride=2

8.3 USING SPECIAL FORTRAN ARRAYS INTRINSIC FUNCTIONS

This type of function will return a value that is dependent on the properties of the object being evaluated.

Chapter 8 Additional Features of Arrays

LBOUND(array)

- Returns ALL of the lower bounds of

array

.

INTEGER :: array ( 1:5, 3:6) WRITE(*,*) "LBOUND is ", LBOUND(array) LBOUND is 1 3

LBOUND(array, dim)

- Returns the lower bound of

dim

of

array

. There is a similar function called

UBOUND .

SHAPE(array)

- Returns the shape of

array

. WRITE(*,*) "SHAPE is ", SHAPE(array) SHAPE is 2 5 4

SIZE(array)

- Returns the number of elements in the array.

WRITE(*,*) "SIZE is ", SIZE(array) SIZE is 20

Example

:

PROGRAM check_array

(Fig.8-7),

FORTRAN provides intrinsic functions for vectors (rank=1 arrays) and matrices (rank=2 arrays)

Table 8.2 (contains the full list; Optional: mask is a LOGICAL array, with element values: TRUE/FALSE)



matmul(matrix_A, matrix_B) dot_product(vector1, vector2) transpose(matrix) maxval(array , mask ) maxloc(array , mask ) minval(array , mask ) minloc(array , mask ) product(array , mask ) sum(array , mask )

From LINEAR ALGEBRA: MATMUL: e.g. A 3x4 x B 4x4 = C 3x4 = 3x where  4

k

  1 DOT_PRODUCT:

b k j

e.g. Vectors

C C C

3 1 . .

. .

. .

X y

X 3x1 = , y 3x1

X

=

y

2

y

3 dot_product computes

i

  1 3

y i

TRANSPOSE: e.g. D 4X3 = A T 3X4 =

a a a a a a a a

MAXVAL:

x

e.g. For array X =

x x

Computes the maximum of x(1) , x(2) , x(3)

a a a a

MINVAL: computes the minimum of x(1), x(2) , x(3) PRODUCT: e.g. for X =

x x x

computes x(1)*x(2)*x(3)

SUM: computes the sum of array elements: x(1) + x(2) +x(3)

An example of matrix and vector multiplication

PROGRAM vectors_and_matrices IMPLICIT NONE INTEGER, DIMENSION(2,3) :: matrix_a INTEGER, DIMENSION(3,2) :: matrix_b INTEGER, DIMENSION(2,2) :: matrix_ab INTEGER, DIMENSION(2) :: vector_c = (/ 1 , 2 /) INTEGER, DIMENSION(3) :: vector_bc ! If all a(i,j) = 1.0, write a = 1.0

! Initialize matrix_a matrix_a(1,1) = 1 matrix_a(1,2) = 2 ! Matrix_a is the matrix:

matrix_a(1,3) = 3 ! [1 2 3] matrix_a(2,1) = 2 ! [2 3 4] matrix_a(2,2) = 3 matrix_a(2,3) = 4 ! Set matrix_b as the transpose of matrix_a matrix_b = TRANSPOSE(matrix_a) !matrix_b is now the matrix [ 1 2 ] ! [ 2 3 ] ! [ 3 4 ] ! Calculate matrix products matrix_ab = MATMUL(matrix_a, matrix_b) ! matrix ab is now the matrix [ 14 20 ] ! [ 20 29 ]

vector_bc = MATMUL(matrix_b, vector_c) !vector_bc is now the vector [5 8 11] END PROGRAM vectors_and_matrices

NOTE 1) The matrix/vectors intrinsic functions help us program very quickly the corresponding linear algebra computations.

e.g. if the program above was programmed with DOloops :

matrix_b = TRANSPOSE(matrix_a) !! Same as the Doloop: DO i = 1,3 DO j = 1,2 matrix_b(i,j) = matrix_a(j,i) ENDDO ENDDO

matrix_ab = MATMUL(matrix_a ,matrix_b) !! Same as the nested Doloop DO i = 1,2 DO j = 1,2 matrix_ab(i ,j) = 0.0

DO k = 1 ,3 matrix_ab(i,j) = matrix_ab(i,j) + matrix_a(i,k) *matrix_b(k,j) ENDDO ENDDO ENDDO

Example program matrix

PROGRAM vectors_and_matrices !! Using array constructors IMPLICIT NONE INTEGER i,j INTEGER , DIMENSION(2 , 3) :: matrix_a = & RESHAPE((/1,2,2,3,3,4/), (/ 2,3 /)) INTEGER, DIMENSION(3,2) :: matrix_b INTEGER, DIMENSION(2,2) :: matrix_ab INTEGER, DIMENSION(2) :: vector_c = (/ 1,2 /) INTEGER, DIMENSION(3) :: vector_bc ! Set matrix_b as the transpose of matrix_a matrix_b = TRANSPOSE(matrix_a) ! Matrix_b is now the matrix [ 1 2 ]

Example program matrix

matrix_ab = MATMUL(matrix_a , matrix_b) !matrix_ab is now the matrix [ 14 20 ] ! [ 20 29 ] vector_bc = MATMUL(matrix_b , vector_c) vector_bc = MATMUL(matrix_b , vector_c) ! Vector_bc is !now the vector : [ 5 8 11 ] WRITE(*,*) 'matrix a = ' WRITE(*, '(3I5) ') ( ( matrix_a(i,j) , j=1,3), i=1,2) WRITE(*, *) 'matrix b =' WRITE(*, '(2I5)' )( ( matrix_b(i,j) , j=1,2) , i=1,3) WRITE(*, *) 'vector b times c' WRITE(*, '(I5)‘) (vector_bc(i) , i=1,3)

Chapter 8 Additional Features of Arrays

Examples of programs: (1) Ex5-8

,

(2) generate (Fig 8-5),

input file: GENDAT

(3)check_array (Fig.8-7).

(4) matrix.f90 (vectors_and_matrices) Note

:

Skip subsections 8.4, 8.5

8.6 ALLOCATABLE ARRAYS

: The arrays studied so far are based on

‘static memory allocation’

(i.e. memory is allocated at the start of the program and remains allocated till the program ends).

Allocatable

arrays are based on

‘dynamic memory allocation’

and allow user to control

allocation / de-allocation

memory for array elements.

of Steps: – Firstly, the allocatable array is specified in a type declaration statement.

– Secondly, space is dynamically allocated for its elements in a separate

allocation statement

, after which the array may be used in the normal way.

– Finally, after the array has been used and is no longer required, the space for the elements is deallocated by a

deallocation statement

.

-Declaration : e.g.

REAL,ALLOCATABLE :: arr1( :,: ) , arr2(: , : ) ALLOCATE(arr1(100,0:10) , arr2(10:30 , 0:10), STAT =status_variable)

NOTE: -After ‘ALLOCATE’ arrays can be used like the standard arrays -Use : STAT = status_variable to check successful allocation.

EXAMPLE PROGRAMS (1)

test_all_array

e.g.

An example of allocation of allocatable arrays, with error checking . … INTEGER :: m=1000,n=100 INTEGER :: error REAL , ALLOCATABLE , DIMENSION(: ,:) :: p INTEGER , ALLOCATABLE , DIMENSION( : ,:) :: q . ALLOCATE(p(n,m), q(10*n,m) , STAT = error) IF(error /= 0) THEN ! Space for p, q could not be allocated WRITE(* ,*) “ Program could not allocate space for p ,q and r ” .

.

STOP END IF ! Space for p , q successfully allocated

-DEALLOCATE (list of allocated arrays, STAT = st_var releases the memory by removing allocated arrays.

e.g.

REAL, ALLOCATABLE, DIMENSION( :) :: varying_array .

.

INTEGER :: i ,n, alloc_error , dealloc_error READ ( * , *) n ! Read maximum size needed DO i = 1 , n ALLOCATE(varying_array( - i : i ) , STAT = alloc_error) IF(alloc_error / = 0) THEN WRITE (* , *) “ Insufficient space to allocate array & &when i = “ , i STOP END IF ! Calculate using varying_array

.

.

.

DEALLOCATE(varying_array, STAT = dealloc_error) IF(dealloc_error / = 0) THEN WRITE (* , *) “ Unexpected deallocation error ” STOP .

.

END IF END DO

- Allocated arrays allow users to have control on the amount of memory space used. This is especially needed for very large n (e.g. if

n

=

10000

then

n 2 = 100000000

) e.g.

PROGRAM space IMPLICIT NONE INTEGER :: n =10000 REAL , ALLOCATABLE , DIMENSION ( : , :) :: a , b ALLOCATE(a(2*n , 6*n)) ! Allocate space for a ! Calculate using a ! Uses 12n 2 elem.’s space .

.

. ! i.e. it needs 12n 2 memory locations for real numbers

.

DEALLOCATE(a) !Free space used by a ALLOCATE(b(3*n , 4*n)) ! Allocate space for b ! Calculate using b ! Uses 12n 2 elem’s space .

.

.

DEALLOCATE (b) !Free space used by b END PROGRAM space Total memory Space = 12n 2 elements space

The same subroutine with standard arrays uses 24

n 2 PROGRAM space IMPLICIT NONE

elements memory space:

INTEGER :: n REAL , DIMENSION(2*n , 6*n) :: a !12n

2 REAL , DIMENSION(3*n , 4*n) :: b !12n

2 elem.’s space elem.’s space !Calculate using a .

.

.

!Calculate using b .

END PROGRAM space ! Total = 24n 2 elem.’ space

NOTE: : There are intrinsic functions which give information on ALLOCATABLE arrays e.g. ‘ALLOCATED’ intrinsic function checks if it was already allocated e.g.

.

REAL , ALLOCATABLE , DIMENSION (:) :: input_data ALLOCATE (input_data (1 : 1000000) , STAT = alloc_stat) …………… IF (ALLOCATED (input_data) ) THEN READ(*,*) input_data DEALLOCATE (input_data) ELSE WRITE (*,*) “warning : array not allocated” END IF Examples: testAllocate1, testAllocate2, testAllocate3

Chapter 9 Additional Features of Procedures

Contents to be covered

:

Section 9.1.1 Explicit-Shape Dummy Arrays Section 9.1.2 Assumed-Shape Dummy Arrays Section 9.2 The SAVE Attribute and Statement Section 9.3 Allocatable Arrays in Procedures Section 9.4 Automatic Arrays in Procedures Remaining Sections are to be Skipped (also skipped) Section 9.1.3 Example 9-1, Fig. 9-1, Fig. 9-2

Chapter 9

9.1 Passing Multidimensional Array to Procedures

The passing of multidimensional (i.e.

rank >=2

) arrays is similar to passing single dimension (i.e.

rank=1

) arrays. However, the procedure needs to know: (1) the number of dimensions (i.e. rank) and (2) the extent of each dimension to use the array. The possible ways to pass arrays are:

Explicit-Shape , Assumed Shape arrays

Chapter 9

9.1.1 Explicit-Shape Arrays

Using this method we pass the array and the extent of each dimension of the array to the subroutine.The extent values are used to declare the size of the array in the procedure (

test example program array2

(

Fig7-6, in Chapter 7

) . Also e.g.

SUBROUTINE process1 (data1, data2, n,m) IMPLICIT NONE INTEGER, INTENT(IN) :: n,m REAL, INTENT(IN), DIMENSION(n,m) :: data1 ! explicit shape REAL, INTENT(OUT), DIMENSION(n,m) :: data2 ! explicit-shape data2= 3.*data1 END SUBROUTINE process1

Chapter 9

e.g.:

Two ways for passing arguments

lower, upper

for explicit-shape arrays:

(1)

As input arguments:

SUBROUTINE explicit (a, b, lower, upper) IMPLICIT NONE INTEGER, INTENT(IN) :: lower, upper REAL, DIMENSION(lower:upper), & INTENT(IN) :: a, b ….

Chapter 9

( 2)

OR, by using a module :

MODULE dataset IMPLICIT NONE INTEGER:: lower, upper lower=10 upper=40 END

MODULE dataset

a, b SUBROUTINE explicit_2(a, b) USE dataset

!

This module includes the lower

!

and upper bounds of a large group

!

of arrays, including a and b

!

These are called lower and upper IMPLICIT NONE REAL, DIMENSION(lower:upper), INTENT(IN) :: & …..

Chapter 9

9.1.2 Assumed-Shape Dummy Arrays

Assumed-shape arrays are declared using colons as placeholders for each subscript of the array. Subroutine must have explicit interface (e.g. must be placed inside a module. Examples in Book (

Skip Example 9-1

),

Example 9-2: PROGRAM assumed_shape (Fig. 9-3).

Chapter 9

Assumed-Shape Dummy Arrays example (not in Book): MODULE test_module CONTAINS SUBROUTINE process1 (data1, data2) REAL, INTENT(IN), DIMENSION(:,:) :: data1 ! Assumed- shape REAL, INTENT(OUT), DIMENSION(:,:) :: data2 ! Assumed-shape data2= 3.*data1 END SUBROUTINE process1 END MODULE test_module

Chapter 9

9.2 The SAVE Attribute and Statement

When the execution of a procedure terminates, the

local variables and local arrays

become undefined (i.e. their values are erased from the memory). But, if the procedure is called many times it is useful to

save

these

values

for the next

CALL

, by using the

SAVE

keyword.

e.g. (1): REAL , SAVE :: sum_x !The value of var sum_x is saved

Also, variables that are initialized are saved even if

SAVE

is not written

e. g. (2): REAL :: sum_x2 = 0. ! The value of var sum_x2 is saved Example (Fig. 9-4, 9-5)

: r

unning_average

,

test_running_average

Chapter 9

9.3 Allocatable Arrays in Procedures (local arrays):

The

allocatable arrays

inside

procedures

are declared the same way as we do in main programs,

but as local arrays

. If we use the

SAVE,

and the procedure is called many times the array-contents will be preserved.

9.5 Allocatable Arrays in Procedures (dummy argument arrays):

Note: We need explicit interface.

1. INTENT(IN) : the array can’t be changed/reallocated/deallocated in the subroutine.

2. INTENT(INOUT) : then the status (‘allocated’ or not) and the data is passed to the subroutine (when it is called). The array can

Chapter 9 be changed/deallocated/re-allocated and the dummy arg. (data) will be passed back to the calling program.

3. INTENT(OUT) : The actual arg. (in calling progr.) will be deallocated (when the subroutine is called) and the subroutine can change/deallocate/re-allocate and the new array will be passed back to the calling program through the arg.

Example: program allocatable_arguments (section 9.5) 9.4 Automatic Arrays in Procedures Automatic Arrays

are

local arrays

inside procedures. They are explicit-shape arrays with non-constant lower/upper bounds.

e.g.

Chapter 9

SUBROUTINE process1 (data1, data2, n,m) IMPLICIT NONE INTEGER, INTENT(IN) :: n,m REAL, INTENT(IN), DIMENSION(n,m) :: data1 ! explicit-shape REAL, INTENT(OUT), DIMENSION(n,m) :: data2 ! explicit-shape REAL, DIMENSION(n,m) :: temp ! Automatic temp=SIN(data1) data2= 3.*data1 + temp END SUBROUTINE process1

Chapter 9 9.4.1 Automatic versus Allocatable arrays: 1) Automatic array is created on entry and removed on exit to/from a procedure (subroutine or function).

2) Allocatable arrays must be allocated/deallocated manually. They can be created/destroyed in (separate procedures). They can be re-sized during calculation.

Note:

Unlike automatic arrays, explicit-shape and assumed-shape arrays can be dummy arg’s and are given

fixed extents

in the calling program.

Chapter 11 Additional intrinsic data types

Objectives:

Understand the different KINDS of a given data type. How to select a specific kind for REAL , INTEGER data.

How to select the precision and range of a real variable in a computer independent manner.

RECALL: REAL Numbers on the computer

REAL

numbers: consider any non-zero real number as a fraction lying between 0.1 and 1.0, called the

mantissa (or fraction)

, which is multiplied or divided by 10 a certain number of times, where this number is called the

exponent

.

a

 10

b

exponent mantissa NOTE: The

mantissa

and the

exponent

are integers numbers.

RECALL: Default (IEEE standard) INTEGER and REAL numbers on a computer The

default

(IEEE standard) INTEGER and REAL numbers on a computer is based on 32-bit representation of numbers. So when we declare data as type

INTEGER

or

REAL

we can only the following numbers

(1)

and

(2)

represented on a computer. The remaining numbers are either

zero

or

+

i

nfinity

or

-

i

nfinity

. For

REAL,

the representation is called

single precision

.

(1)

An

integer

is a whole number, is always stored

exactly

in the computer’s memory, and has a (relatively) limited range : (between about -2*10 9 and +2*10 9 on a typical 32-bit computer).

(2)

A

real

number, on the other hand, is stored as a floating-point number, is stored as an

approximation

to a fixed number of significant digits and it has a very large range : (typically in the range between about -10 38 and +10 38 and

accurate

in the first

seven

or

eight

significant digits on the same 32-bit computer).

Chapter 11 Additional intrinsic data types

11.1 Alternative kinds of REAL data type

: In Fortran, there is a

double precision

(

REAL

data type ), which is based on a representation on the computer with 64-bits (i.e. 8 bytes). For this, the mantissa (or fraction of the number) have 14 (depending on the computer maybe: 15 or 16) accurate decimal digits and the exponent range is [-308 , 308] (for base 10). The content of this chapter is (in part) devoted in explaining how to use

double precision

in

Fortran 90/95 programs

.

Chapter 11 Additional intrinsic data types

11.1.1 Kinds of REAL constants and variables:

The KIND parameter can be used in the REAL declaration to specify single (or double) precision. e.g. REAL(KIND=1) :: value_1 REAL(KIND= 4) :: value_2 REAL(KIND=8 ), DIMENSION(20) :: array REAL(4) :: temp ! better to use ‘KIND=‘ We often call the variables declared (with KIND)

parametrized variables

. In general, KIND=1 or 4 is single and KIND=2 or 8 is double precision. But, it all depends on the computer/compiler.

Chapter 11 Additional intrinsic data types In order, to make the Fortran programs

portable

among computers we use a

named constant

(instead of 4 or 8) in KIND: e.g.

INTEGER, PARAMETER :: SINGLE = 4 INTEGER, PARAMETER :: DOUBLE = 8 REAL(KIND= SINGLE) :: value_1 REAL(KIND= DOUBLE), DIMENSION(2) :: array REAL(SINGLE) :: temp

We can also declare constants: e.g.

34.

! default kind,

single precision 34._4

34.E3

! Valid if KIND=4 is valid; ! default kind,

please avoid this way single precision

!

1234.56789_DOUBLE

! Valid, if DOUBLE is declared, as above

Chapter 11 Additional intrinsic data types 3.0E0 ! default kind,

single precision

constant 3.0D0 !

double precision

constant We can find out the kind of precision on a computer by using the intrinsic function KIND. e.g.

PROGRAM kinds (Figure 11-1).

Chapter 11 Additional intrinsic data types

11.1.3 Selecting a precision (in a computer/processor independent manner)

We use an

intrinsic function

SELECTED_REAL_KIND : kind_number = SELECTED_REAL_KIND (p=precision, & r=range) ! precision is the no. of decimal digits of precision, i.e. accuracy ! range is the exponents range in powers of 10.

! e.g. kind_number = SELECTED_REAL_KIND (p=6, r=37) This function gives: -1 (

precision is not available

), -2 (

range is not available

), -3 (

precision and range are not available

).

Note

: what is available depends on the computer/compiler.

Chapter 11 Additional intrinsic data types

Table 11-2

: Kind related functions return the following:

SELECTED_REAL_KIND (p, r)

! the smallest kind of REAL with p ! decimal digits of precision and

range

: [-10^r, 10^r]

SELECTED_INT_KIND (r)

! Smallest kind of INTEGER ! with

range:

[-10^r, 10^r]

KIND(X)

! the KIND for the var. or const. X

PRECISION(X)

! the decimal digits of precision for the var. or const. X

RANGE(X)

! the decimal exponent range for the var. or const. X

Example: PROGRAM selected_kinds

(Figure 11-2)

Chapter 11 Additional intrinsic data types

11.1.4 Mixed-precision (single/double precision) arithmetic Example (operation + in single/double precision): 1./3 + 1./3

1.

1.D0/3. + 1/3 = 3.333333333333333 E-001 ! 1/3 =0, here 2.

1./3. + 1.D0/3. = 6.666666633333333 E-001 ! 1./3., is single precision 3.

1.D0/3. + 1./3.D0 = 6.666666666666666 E-001

Example

: Constants must be declared/written in double precision to actually have a double precision accuracy. e.g.

PROGRAM test_initial INTEGER, PARAMETER :: DBL =SELECTED_REAL_KIND(p=13) REAL(KIND=DBL) :: a1= 6.666666666666666

Chapter 11 Additional intrinsic data types

REAL(KIND=DBL) :: a2= 6.666666666666666_DBL

WRITE(*,*) a1,a2 END PROGRAM test_initial

The program gives the results for a1 and a2: 6.666666507720947 6.666666666666666

(i.e. a1 is in single-precision , because it was not declared/written as double!).

11.1.5 Higher precision intrinsic functions

:

y=DBLE(x)

! converts the x to double precision y. But

DBLE(a1)

will not ! Add lost digits if constant

a1

did not have

_DBL,

like

a2 Note: y=REAL(x)

! converts the x to single precision y.

Chapter 11 Additional intrinsic data types

11.1.6 When should one use high precision (i.e. double precision) instead of standard REAL (?)

e.g. REAL :: x, y, a1, a2

1. When the values are : x < -10^(38) or x > 10^(38) 2. When we have: x+y (or x-y) with |x| much smaller than |y| e.g. x = 3.25 and y= 1000000.0 then: Single-precision: x+y= 1000003.0 Double-precision: x+y= 1000003.25

3. When the |a1| and |a2| are about the same and we compute (a1-a2) .

Chapter 11 Additional intrinsic data types

Example : Numerical calculation of derivatives.

e.g. PROGRAM diff (Figure 11-4). Calculate the derivative of f(x) = 1/x , at x= 0.15 .

Note: (1) the true derivative of f(x) is: -1/(x*x) = 44.444444444 .. ; (2) using the approximation formula: (f(x+d_x) –f(x) )/d_x Note

: The program uses

d_x

= 1.000E-1 ,.., 1.000E-10 . When d_x is smaller than 1.000E-5 both single and double precision solutions are accurate. But for d_x > 1.000E-5 only the double precision solution is accurate .

Chapter 11 Additional intrinsic data types

SKIP : section 11.1.7

11.2 Alternate lengths of the INTEGER data type

The intrinsic function

kind_number = SELECTED_INT_KIND(e)

(where range is the exponent e for the range [-10^e , 10^e] , e.g. e=3,9,12,30) gives the

integer parameter

, to be used in a declaration for INTEGER type.

e.g. the program:

PROGRAM test_selected_int_kinds Note

:

e.g e=9 is the default range

the range is

+ infinity

or for

– infinity

.

INTEGER

on Intel Pentiums : [-2147483648 , 2147483648 ]. A number outside

Chapter 11 Additional intrinsic data types Example of declaration. INTEGER, PARAMETER :: SHORT= SELECTED_INT_KIND(3) INTEGER, PARAMETER :: LONG = SELECTED_INT_KIND(9) INTEGER (KIND=SHORT) :: i1 INTEGER (KIND=LONG) :: i2 INTEGER :: 34, 34_4, 24_LONG ! the constant 34_4 is valid if 4 is a legal kind of integer

SKIP: the remaining sections