Docstoc

Fortran 90 Tutorial_4

Document Sample
Fortran 90 Tutorial_4 Powered By Docstoc
					FUNCTIONS AND MODULES

SELECT THE TOPICS YOU WISH TO REVIEW:
  Functions and Modules

         Designing Functions

         Functions Examples

         Common Problems

         Using Functions

         Argument Association

         Where Do My Functions Go?

         Scope Rules      Programming Examples:

         Computing Cubes

         Computing Means - Revisited

         Cm and Inch Conversion

         Heron'a Formula for Computing Triangle Area - Revisited

         Computing the Combinatorial Coefficient

         Computing Square Roots with Newton's Method

         Designing a Greatest Common Divisor Function

         Finding All Prime Numbers in the Range of 2 and N - Revisited

         The Bisection (Bozano) Equation Solver

         A Brief Introduction to Modules

         What is a Module?

         How to Use a Module?

         Compile Programs with Modules

         Programming Example 1: Factorial and Combinatorial Coefficient

         Programming Example 2: Trigonometric Functions Using Degree
          A Little Privacy - PUBLIC/PRIVATE

           External Functions and Interface Blocks

          Interface Blocks

          Programming Example 1: Cm and Inch Conversion

          Programming Example 2: Heron'a Formula for Computing Triangle Area

           Download my course overheads


DESIGNING FUNCTIONS

SYNTAX
In addition to intrinsic functions, Fortran allows you to design your own functions. A Fortran function, or more
precisely, a Fortran function subprogram, has the following syntax:

type FUNCTION function-name (arg1, arg2, ..., argn)
   IMPLICIT NONE
   [specification part]
   [execution part]
   [subprogram part]
END FUNCTION function-name

Here are some elaborations of the above syntax:

        The first line of a function starts with the keyword FUNCTION. Before FUNCTION, the type gives the type
         of the function value (i.e., INTEGER, REAL, LOGICAL and CHARACTER) and after FUNCTION is the name
         you assign to that function.
        Following the function-name, there is a pair of parenthesis in which a number of arguments arg1, arg2, ...,
         argn are separated with commas. These arguments are referred to as formal arguments. Formal
         arguments must be variable names and cannot be expressions. Here are a examples:
              1. The following is a function called Factorial. It takes only one formal argument n and returns an
                   INTEGER as its function value.
                        INTEGER FUNCTION Factorial(n)
              1. The following is a function called TestSomething. It takes three formal arguments a, b and c, and
                   returns a LOGICAL value (i.e., .TRUE. or .FALSE.) as its function value.
                        LOGICAL FUNCTION TestSomething(a, b, c)
        A function must be ended with END FUNCTION followed by the name of that function.
        Between FUNCTION and END FUNCTION, there are the IMPLICIT NONE, specification part, execution part
         and subprogram part. These are exactly identical to that of a PROGRAM.

If a function does not need any formal argument, it can be written as
type FUNCTION function-name ()
      IMPLICIT NONE
      [specification part]
      [execution part]
      [subprogram part]
END FUNCTION function-name
where arg1, arg2, ..., argn are left out. But, the pait of parenthesis must be there.




SEMANTICS
The meaning of a function is very simple:

        A function is a self-contained unit that receives some "input" from the outside world via its formal
         arguments, does some computations, and then returns the result with the name of the function.
        Thus, since the function returns its result via the name of the function, somewhere in the function there
         must exist one or more assignment statements like the following:
           function-name = expression

         where the result of expression is stored to the name of the function.

         However, the name of the function cannot appear in the right-hand side of any
         expression. That is, the name of the function, used as a variable name, can only appear in
         the left-hand side of an expression. This is an artificial restriction in this course only.

        A function receives its input values from formal arguments, does computations, and saves the result in its
         name. When the control of execution reaches END FUNCTION, the value stored in the name of the
         function is returned as the function value.
        To tell the function about the types of its formal arguments, all arguments must be declared with a new
         attribute INTENT(IN). The meaning of INTENT(IN) indicates that the function will only take the value from
         the formal argument and must not change its content.
        Any statements that can be used in PROGRAM can also be used in a FUNCTION.




FUNCTIONS EXAMPLES


Here are a few examples of functions:

        The following function has a name Sum and three formal arguments a, b and c. It returns an INTEGER
         function value. The INTEGER, INTENT(IN) part indicates that the function takes its input value from its
         three formal argument. Then, the function uses the value of these formal arguments to compute the sum
         and stores in Sum, the name of the function. Since the next statement is END FUNCTION, the function
         returns the value stored in Sum.

           INTEGER FUNCTION Sum(a, b, c)
              IMPLICIT NONE

                INTEGER, INTENT(IN) :: a, b, c

              Sum = a + b + c
           END FUNCTION Sum
    If the value supplied to a, b and c are 3, 5, and -2, respectively, Sum will receive 6 (=3+5+(-2)) and the
    function returns 6.

   The following function has a name Positive with a REAL formal argument. If the argument is positive, the
    function returns .TRUE.; otherwise, the function returns .FALSE.

      LOGICAL FUNCTION Positive(a)
         IMPLICIT NONE

          REAL, INTENT(IN) :: a

         IF (a > 0.0) THEN
            Positive = .TRUE.
         ELSE
            Positive = .FALSE.
         END IF
      END FUNCTION Positive

    The above function can be made much shorter by using LOGICAL assignment. In the following, if a > 0.0 is
    true, .TRUE. is stored to Positive; otherwise, Positive receives .FALSE.

    LOGICAL FUNCTION Positive(a)
       IMPLICIT NONE

        REAL, INTENT(IN) :: a

       Positive = a > 0.0
    END FUNCTION Positive

   The following function, LargerRoot, takes three REAL formal arguments and returns a REAL function
                                                                2
    value. It returns the larger root of a quadratic equation ax + bx + c = 0.

      REAL FUNCTION LargerRoot(a, b, c)
         IMPLICIT NONE

          REAL, INTENT(IN) :: a
          REAL, INTENT(IN) :: b
          REAL, INTENT(IN) :: c
          REAL             :: d, r1, r2

          d = SQRT(b*b - 4.0*a*c)
          r1= (-b + d) / (2.0*a)
          r2= (-b - d) / (2.0*a)
          IF(r1 >= r2) THEN
            LargerRoot = r1
         ELSE
            LargerRoot = r2
         END IF
      END FUNCTION LargerRoot

    The above example shows that you can declare other variables such as d, r1 and r2 if they are needed.
   The following function, Factorial(), has only one INTEGER formal argument n >= 0, and computes and
    returns the factorial of n, n!.

      INTEGER FUNCTION Factorial(n)
         IMPLICIT NONE

          INTEGER, INTENT(IN) :: n
          INTEGER             :: i, Ans

         Ans = 1
         DO i = 1, n
            Ans = Ans * i
         END DO
         Factorial = Ans
      END FUNCTION

    Note that the function name Factorial is not used in any computation. Instead, a new INTEGER variable is
    used for computing n!. The final value of Ans is stored to Factorial before leaving the function.

    If Factorial is involved in computation like the following:

    INTEGER FUNCTION Factorial(n)
       IMPLICIT NONE

        INTEGER, INTENT(IN) :: n
        INTEGER             :: i

       Factorial = 1
       DO i = 1, n
          Factorial = Factorial * i
       END DO
    END FUNCTION

    then there is a mistake although the above program looks normal. The reason is that the function name
    cannot appear in the right-hand side in any expression of that function.

   The following function GetNumber() does not have any formal arguments and returns an INTEGER
    function value. This function has a DO-loop which keeps asking the user to input a positive number. The
    input is read into the function name. If this value is positive, then EXIT and the function returns the value
    in GetNumber. Otherwise, the loop goes back and asks the user again for a new input value.

      REAL FUNCTION GetNumber()
         IMPLICIT NONE

          DO
            WRITE(*,*) 'A positive real number --> '
            READ(*,*)   GetNumber
            IF (GetNumber > 0.0) EXIT
            WRITE(*,*) 'ERROR. Please try again.'
         END DO
         WRITE(*,*)
      END FUNCTION GetNumber
COMMON PROBLEMS

        Forget the type of a FUNCTION.

        FUNCTION DoSomething(a, b)
           IMPLICIT NONE

            INTEGER, INTENT(IN) :: a, b

           DoSomthing = SQRT(a*a + b*b)
        END FUNCTION DoSomthing

      If there is no type, you will not be able to determine if the returned value is an INTEGER, a REAL or
      something else.

        Forget INTENT(IN).

        REAL FUNCTION DoSomething(a, b)
           IMPLICIT NONE

            INTEGER :: a, b

           DoSomthing = SQRT(a*a + b*b)
        END FUNCTION DoSomthing

      Actually, this is not an error. But, without INTENT(IN), our Fortran compiler will not be able to check many
      potential errors for you.

     Change the value of a formal argument declared with INTENT(IN).

        REAL FUNCTION DoSomething(a, b)
           IMPLICIT NONE

            INTEGER, INTENT(IN) :: a, b

           IF (a > b)        THEN
              a = a -        b
           ELSE
              a = a +        b
           END IF
           DoSomthing        = SQRT(a*a + b*b)
        END FUNCTION         DoSomthing

      Since formal argument a is declared with INTENT(IN), its value cannot be changed.

     Forget to store a value to the function name.

        REAL FUNCTION DoSomething(a, b)
           IMPLICIT NONE
             INTEGER, INTENT(IN) :: a, b
             INTEGER             :: c

            c = SQRT(a*a + b*b)
         END FUNCTION DoSomthing

       When the execution of this function reaches END FUNCTION, it returns the value stored in DoSomething.
       However, in the above, since there is no value ever stored in DoSmething, the returned value could be
       anything (i.e., a garbage value).

   Function name is used in the right-hand side of an expression.

   REAL FUNCTION DoSomething(a, b)
      IMPLICIT NONE

            INTEGER, INTENT(IN) :: a, b

           DoSomething = a*a + b*b
           DoSomething = SQRT(DoSomething)
        END FUNCTION DoSomthing

       In the above, function name DoSomething appears in the right-hand side of the second expression. This is
       a mistake. Only a special type of functions, RECURSIVE functions, could have their names in the right-hand
       side of expressions.

      Only the most recent value stored in the function name will be returned..

         REAL FUNCTION DoSomething(a, b)
            IMPLICIT NONE

             INTEGER, INTENT(IN) :: a, b

            DoSomething = a*a + b*b
            DoSomething = SQRT(a*a - b*b)
         END FUNCTION DoSomthing

       In the above, the value of SQRT(a*a-b*b) rather than the value of a*a + b*b is returned. In fact, this is
       obvious. Since the name of a function can be considered as a special variable name, the second
       assignment overwrites the previous value.


USING FUNCTIONS


The way of using a user-defined function is exactly identical to that of using Fortran intrinsic
functions. One can use a function in an expression and in a WRITE. Suppose we have the
following function:

REAL FUNCTION Average(x, y, z)
   IMPLICIT NONE
    REAL, INTENT(IN) :: x, y, z

     Average = (x + y + z) / 3.0
END FUNCTION Average
This function takes three REAL formal arguments and returns their average.

To use this function, one needs to supply values to the formal arguments. For example, one could
write the following:

..... = ..... + Average(1.0, 2.0, 3.0) * .....
The above expression involves the use of function Average. Since this function has three formal arguments, three
values must be presented to the function. Here, the values are 1.0, 2.0 and 3.0 and the returned value is the
average of these three numbers (i.e., 2.0). The values or expressions used to invoke a function are referred to as
actual arguments.

Please keep the following important rules in mind:

    The number of formal arguments and actual arguments must be equal.

        ..... = ... + Average(1.0, 2.0, 3.0, 4.0) * ....
        ..... = ... - Average(3.0, 6.0) + .....

        The first line has more actual arguments than the number of formal arguments, and the second line has
        less actual arguments than the number of formal arguments.

    The type of the corresponding actual and formal arguments must be identical.

    WRITE(*,*) Average(1, 2.5, 3.7)

        In the above example, the first actual argument is an INTEGER which doe not match with the type of the
        first formal argument. Thus, it is incorrect.

    The actual arguments can be constants, variables and even expressions.

    REAL :: a = 1.0, b = 2.0, c = 3.0

    ..... = ... + Average(1.0, 2.0, 3.0) + .....
    ..... = ... + Average(a, b, c) + .....
    ..... = ... + Average(a+b, b*c, (b+c)/a) + .....

        In the above, the first line shows the use of constants as actual arguments. The second line uses variables,
        while the third uses expression. In the third line, the result of evaluating a+b, b*c and (b+c)/a will be
        supplied to the three formal arguments of function Average().

Please continue with the next page on argument association
ARGUMENT ASSOCIATION


When using a function, the values of actual arguments are passed to their corresponding formal
arguments. There are some important rules:

       If an actual argument is an expression, it is evaluated and the result is saved into a temporary location.
        Then, the value in this temporary location is passed.
       If an actual argument is a constant, it is considered as an expression. Therefore, its value is saved to a
        temporary location and then passed.
       If an actual argument is a variable, its value is taken and passed to the corresponding formal argument.
       If an actual argument is a variable enclosed in a pair of parenthesis like (A), then this is an expression and
        its value is evaluated and saved to a temporary location. Then, this value (in the temporary location) is
        passed.
       For a formal argument declared with INTENT(IN), any attempt to change its value in the function will
        cause a compiler error. The meaning of INTENT(IN) is that this formal argument only receives a value
        from its corresponding actual argument and its value cannot be changed.
       Based on the previous point, if a formal argument is declared without using INTENT(IN), its value can be
        changed. For writing functions, this is not a good practice and consequently all formal arguments should
        be declared with INTENT(IN).



EXAMPLE
We shall use the following function to illustrate the above rules. Function Small() takes three formal arguments x,
y and z and returns the value of the smallest.

INTEGER :: a, b, c                           INTEGER FUNCTION Small(x, y, z)
                                                IMPLICIT NONE
a = 10                                          INTEGER, INTENT(IN) :: x, y, z
b = 5
c = 13                                          IF (x <= y        .AND. x <= z) THEN
WRITE(*,*)      Small(a,b,c)                       Small =        x
WRITE(*,*)      Small(a+b,b+c,c)                ELSE IF (y        <= x .AND. y <= z) THEN
WRITE(*,*)      Small(1, 5, 3)                     Small =        y
WRITE(*,*)      Small((a),(b),(c))              ELSE
                                                   Small =        z
                                                END IF
                                             END FUNCTION         Small

In the first WRITE, all three arguments are variable and their values are sent to the
corresponding formal argument. Therefore, x, y and z in function Small() receives 10, 5 and 13,
respectively. The following diagram illustrate this association:
In the second WRITE, the first and second arguments are expression a+b, b+c. Thus, they are
evaluated yielding 15 and 18. These two values are saved in two temporary locations and passed
to function Small() along with variable c. Thus, x, y and z receive 15, 18 and 13 respectively.
This is illustrated as follows. In the figure, squares drawn with dotted lines are temperary
locations that are created for passing the values of arguments.




In the third WRITE, since all actual arguments are constants, their values are "evaluated" to
temporary locations and then sent to x, y and z.




In the fourth WRITE, since all three actual arguments are expressions, they are evaluated and
stored in temporary locations. Then, these values are passed to x, y and z.
Thus, x, y and z receive 10, 5 and 13, respectively.

Note that while the formal arguments of function Small() receive the same values using the
first and the fourth lines, the argument associations are totally different. The first line has
the values in variables passed directly and the the fourth line evaluates the expressions into
temporary locations whose values are passed. This way of argument association does not
have impact on functions (since you must use INTENT(IN)), it will play an important role
in subroutine subprograms.

WHERE DO MY FUNCTIONS GO?


Now you have your functions. Where do they go? There are two places for you to add these
functions in. They are either internal or external. This page only describes internal functions.
See external subprogram for the details of using external functions.

Long time ago, we mentioned the structure of a Fortran program. From there, we know that the
last part of a program is subprogram part. This is the place for us to put the functions. Here is the
syntax:

PROGRAM program-name
IMPLICIT NONE
[specification part]
[execution part]
CONTAINS
    [your functions]
END PROGRAM program-name
In the above, following all executable statements, there is the keyword CONTAINS, followed by all of your
functions, followed by END PROGRAM.

From now on, the program is usually referred to as the main program or the main program
unit. A program always starts its execution with the first statement of the main program. When a
function is required, the control of execution is transfered into the corresponding function until
the function completes its task and returns a function values. Then, the main program continues
its execution and uses the returned function value for further computation.


EXAMPLES

       THE FOLLOWING COMPLETE PROGRAM "CONTAINS" ONE FUNCTION, AVERAGE(). THE EXECUTION PART
        CONSISTS OF THREE STATEMENTS, A READ, AN ASSIGNMENT AND A WRITE. THESE THREE
        STATEMENTS MUST BE PLACED BEFORE THE KEYWORD CONTAINS. N OTE ALSO THAT END
        PROGRAM MUST BE THE LAST LINE OF YOUR PROGRAM .

           PROGRAM Avg
              IMPLICIT NONE
              REAL :: a, b, c, Mean
            READ(*,*) a, b, c
            Mean = Average(a, b, c)
            WRITE(*,*) a, b, c, Mean

        CONTAINS
           REAL FUNCTION Average(a, b, c)
              IMPLICIT NONE
              REAL, INTENT(IN) :: a, b, c
              Average = (a + b + c) / 3.0
           END FUNCTION Average
        END PROGRAM Avg

      THE FOLLOWING PROGRAM "CONTAINS" TWO FUNCTIONS LARGE() AND GEOMEAN(). THE ORDER OF
       THESE FUNCTIONS ARE UNIMPORTANT .

        PROGRAM TwoFunctions
           IMPLICIT NONE
           INTEGER :: a, b, BiggerOne
           REAL    :: GeometricMean

            READ(*,*) a, b
            BiggerOne = Large(a,b)
            GeometricMean = GeoMean(a,b)
            WRITE(*,*) 'Input = ', a, b
            WRITE(*,*) 'Larger one = ', BiggerOne
            WRITE(*,*) 'Geometric Mean = ', GeometricMean

        CONTAINS
           INTEGER FUNCTION Large(a, b)
              IMPLICIT NONE
              INTEGER, INTENT(IN) :: a, b
              IF (a >= b) THEN
                 Large = a
              ELSE
                 Large = b
              END IF
           END FUNCTION Large

           REAL FUNCTION GeoMean(a, b)
              IMPLICIT NONE
              INTEGER, INTENT(IN) :: a, b
              GeoMean = SQRT(REAL(a*b))
           END FUNCTION GeoMean
        END PROGRAM TwoFunctions


A IMPORTANT NOTE
                    Although in general a function can "contain" other functions,
                    an internal function CANNOT contain any other functions.
In the following example, the main program "contains" an internal function InternalFunction(),
which in turn contains another internal function Funct(). This is incorrect, since an internal
function cannot contain another internal function. In other words, the internal functions of a main
program must be on the same level.

PROGRAM Wrong
    IMPLICIT NONE
         .........
CONTAINS
    INTEGER FUNCTION InternalFunction(.....)
         IMPLICIT NONE
             ........
    CONTAINS
         REAL FUNCTION Funct(.....)
             IMPLICIT NONE
                 ........
         END FUNCTION Funct
    END FUNCTION InternalFunction
END PROGRAM Wrong
Please continue with the next important topic about scope rules.




SCOPE RULES


Since a main program could contain many functions and in fact a function can contain other
functions (i.e., nested functions), one may ask the following questions:

    1.   Could a function use a variable declared in the main program?
    2.   Could a main program use a variable declared in one of its function?

The scope rules answer these questions. In fact, scope rules tell us if an entity (i.e., variable, parameter and
function) is "visible" or accessible at certain places. Thus, places where an entity can be accessed or visible is
referred to the scope of that entity.

The simplest rule is the following:

                                        The scope of an entity is the program or function
                         Scope Rule 1
                                        in which it is declared.



Therefore, in the following, the scope of parameter PI and variables m and n is the main
program, the scope of formal argument k and REAL variables f and g is function Funct1(), and
the scope of formal arguments u and v is function Funct2().

PROGRAM Scope_1
   IMPLICIT NONE
   REAL, PARAMETER :: PI = 3.1415926
   INTEGER          :: m, n
      ...................
CONTAINS
   INTEGER FUNCTION Funct1(k)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: k
      REAL                  :: f, g
         ..........
   END FUNCTION Funct1

   REAL FUNCTION Funct2(u, v)
      IMPLICIT NONE
      REAL, INTENT(IN) :: u, v
         ..........
   END FUNCTION Funct2
END PROGRAM Scope_1

There is a direct consequence of Scope Rule 1. Since an entity declared in a function has a scope
of that function, this entity cannot be seen from outside of the function. In the above example,
formal argument k and variables f and g are declared within function Funct1(), they are only
"visible" in function Funct1() and are not visible outside of Funct1(). In other words, since k, f
and g are not "visible" from the main program and function Funct2(), they cannot be used in the
main program and function Funct2(). Similarly, the main program and function Funct1() cannot
use the formal arguments u and v and any entity declared in function Funct2().


LOCAL ENTITIES
Due to the above discussion, the entities declared in a function or in the main program are said local to that
function or the main program. Thus, k, f and g are local to function Funct1(), u and v are local to function Funct2(),
and PI, m and n are local to the main program.




GLOBAL E NTITIES
Given a function f(), entities that are declared in all containing functions or the main program are said global to f().
In the above example, since variables m and n are declared in the main program, they are global to Funct1() and
function Funct2(). However, variables f and g are not global to function Funct2(), since Funct1() does not contain
Funct2(). Similarly, formal arguments u and v are not global to function Funct1().

This comes the second scope rule:

                                      A global entity is visible to all contained functions,
                       Scope Rule 2
                                      including the function in which that entity is declared.



Continue with the above example, since m and n are global to both functions Funct1() and
Funct2(), they can be used in these two functions.
PROGRAM Scope_2
   IMPLICIT NONE
   INTEGER :: a = 1, b = 2, c = 3

    WRITE(*,*)        Add(a)
    c = 4
    WRITE(*,*)        Add(a)
    WRITE(*,*)        Mul(b,c)

CONTAINS
   INTEGER FUNCTION Add(q)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: q
      Add = q + c
   END FUNCTION Add

   INTEGER FUNCTION Mul(x, y)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: x, y
      Mul = x * y
   END FUNCTION Mul
END PROGRAM Scope_2

In the above program, variables a, b and c are global to both functions Add() and Mul(). Therefore, since variable c
used in function Add() is global to Add(), expression q + c means computing the sum of the value of the formal
argument q and the value of global variable c. Therefore, the first WRITE produces 4 (= 1 + 3). Before the second
WRITE, the value of c is changed to 4 in the main program. Hence, the second WRITE produces 5 (= 1 + 4). The
third WRITE produces 8 (= 2 * 4).

Thus, the first two WRITEs produce different results even though their actual arguments are the
same! This is usually refereed to as a side effect. Therefore, if it is possible, avoid using global
variables in internal functions.

Let us continue with the above example. To remove side effect, one could add one more
argument to function Add() for passing the value of c.

PROGRAM Scope_2
   IMPLICIT NONE
   INTEGER :: a = 1, b = 2, c = 3

    WRITE(*,*)        Add(a, c)
    c = 4
    WRITE(*,*)        Add(a, c)
    WRITE(*,*)        Mul(b,c)

CONTAINS
   INTEGER FUNCTION Add(q, h)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: q, h
      Add = q + h
   END FUNCTION Add

    INTEGER FUNCTION Mul(x, y)
       IMPLICIT NONE
      INTEGER, INTENT(IN) :: x, y
      Mul = x * y
   END FUNCTION Mul
END PROGRAM Scope_2


WHAT IF THERE ARE NAME CONFLICTS ?
Frequently we may have a local entity whose name is identical to the name of a global entity. To resolve this name
conflict, we need the following new scope rule:

                                   An entity declared in the scope of another entity is always
                    Scope Rule 3
                                   a different entity even if their names are identical.



In the program below, the main program declares a variable i, which is global to function Sum().
However, i is also declared in function Sum(). According to Scope Rule 3 above, these two is are
two different entities. More precisely, when the value of Sum()'s i is changed, this change will
not affect the i in the main program and vice versa. This would save us a lot of time in finding
variables with different names.

PROGRAM Scope_3
   IMPLICIT NONE
   INTEGER :: i, Max = 5

    DO i = 1, Max
       Write(*,*)         Sum(i)
    END DO

CONTAINS

   INTEGER FUNCTION Sum(n)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: n
      INTEGER             :: i, s
      s = 0
      DO i = 1, n
         s = s + i
      END DO
      Sum = s
   END FUNCTION Sum
END PROGRAM Scope_3




Programming Examples:
COMPUTING CUBES

PROBLEM STATEMENT
Write a program to compute the cubes of 1, 2, 3, ..., 10 in both INTEGER and REAL types. It is
required to write a function intCube() for computing the cube of an integer and a function
realCube() for computing the cube of a real.


SOLUTION
!   -----------------------------------------------------
!      This program display the cubes of INTEGERs and
!   REALs. The cubes are computed with two functions:
!   intCube() and realCube().
!   -----------------------------------------------------

PROGRAM Cubes
   IMPLICIT   NONE

     INTEGER, PARAMETER :: Iterations = 10
     INTEGER            :: i
     REAL               :: x

     DO i = 1, Iterations
        x = i
        WRITE(*,*) i, x, intCube(i), realCube(x)
     END DO

CONTAINS

! -----------------------------------------------------
! INTEGER FUNCTION intCube() :
!    This function returns the cube of the argument.
! -----------------------------------------------------

     INTEGER FUNCTION intCube(Number)
        IMPLICIT  NONE

        INTEGER, INTENT(IN) :: Number

        intCube = Number*Number*Number
     END FUNCTION intCube

! -----------------------------------------------------
! REAL FUNCTION realCube() :
!    This function returns the cube of the argument.
! -----------------------------------------------------

     REAL FUNCTION realCube(Number)
        IMPLICIT NONE

        REAL, INTENT(IN) :: Number

        realCube = Number*Number*Number
     END FUNCTION realCube

END PROGRAM Cubes
Click here to download this program.
PROGRAM INPUT AND OUTPUT
The following is the output from the above program.

 1,     1., 1, 1.
 2,     2., 8, 8.
 3,     3., 27, 27.
 4,     4., 64, 64.
 5,     5., 125, 125.
 6,     6., 216, 216.
 7,     7., 343, 343.
 8,     8., 512, 512.
 9,     9., 729, 729.
 10,     10., 1000, 1000.


DISCUSSION

       Functions intCube() and realCube() are identical except for the types of formal arguments and returned
        values.
       The main program has a DO-loop that iterates Iterations times (this is a PARAMETER, or an alias, of 10).
        Variable x holds the real value of integer variable i.


COMPUTING MEANS - REVISITED

PROBLEM STATEMENT

The arithmetic, geometric and harmonic means of three positive numbers are defined by the
following formulas:




Write a program to read three positive numbers and use three functions to compute the
arithmetic, geometric and harmonic means.


SOLUTION
! ----------------------------------------------------------
!    This program contains three functions for computing the
! arithmetic, geometric and harmonic means of three REALs.
! ----------------------------------------------------------
PROGRAM ComputingMeans
   IMPLICIT NONE

     REAL   :: a, b, c

     READ(*,*)     a, b, c
     WRITE(*,*)    'Input: ', a, b, c
     WRITE(*,*)
     WRITE(*,*)   'Arithmetic mean = ', ArithMean(a, b, c)
     WRITE(*,*)   'Geometric mean = ', GeoMean(a, b, c)
     WRITE(*,*)   'Harmonic mean   = ', HarmonMean(a, b, c)

CONTAINS

!   ----------------------------------------------------------
!   REAL FUNCTION ArithMean() :
!      This function computes the arithmetic mean of its
!   three REAL arguments.
!   ----------------------------------------------------------

     REAL FUNCTION ArithMean(a, b, c)
        IMPLICIT NONE

        REAL, INTENT(IN) :: a, b, c

        ArithMean = (a + b + c) /3.0
     END FUNCTION ArithMean

!   ----------------------------------------------------------
!   REAL FUNCTION GeoMean() :
!      This function computes the geometric mean of its
!   three REAL arguments.
!   ----------------------------------------------------------

     REAL FUNCTION GeoMean(a, b, c)
        IMPLICIT NONE

        REAL, INTENT(IN) :: a, b, c

        GeoMean = (a * b * c)**(1.0/3.0)
     END FUNCTION GeoMean

!   ----------------------------------------------------------
!   REAL FUNCTION HarmonMean() :
!      This function computes the harmonic mean of its
!   three REAL arguments.
!   ----------------------------------------------------------

     REAL FUNCTION HarmonMean(a, b, c)
        IMPLICIT NONE

        REAL, INTENT(IN) :: a, b, c

        HarmonMean = 3.0 / (1.0/a + 1.0/b + 1.0/c)
     END FUNCTION HarmonMean

END PROGRAM       ComputingMeans
Click                here                 to                download                   this                program.



PROGRAM INPUT AND OUTPUT
The following is the output from the above program.

Input: 10.,        15.,     18.

Arithmetic mean = 14.333333
Geometric mean = 13.9247675
Harmonic mean   = 13.5


DISCUSSION

        Each of these functions is simple and does not require further explanation.
        Note that the main program and all three functions use the same names a, b and c. By Scope Rule 3, they
         are all different entities. That is, when function ArithMean() is using a, it is using its own local formal
         argument rather than the global variable a declared in the main program.


CM AND INCH CONVERSION

PROBLEM STATEMENT

It is known that 1 cm is equal to 0.3937 inch and 1 inch is equal to 2.54 cm. Write a program to
convert 0, 0.5, 1, 1.5, ..., 8, 8.5, 9, 9.5, and 10 from cm to inch and from inch to cm.


SOLUTION
!   ---------------------------------------------------------------
!      This program "contains" two REAL functions:
!           (1) Cm_to_Inch() takes a real inch unit and converts
!                it to cm unit, and
!           (2) Inch_to_cm() takes a real cm unit and converts it
!                to inch unit.
!   The main program uses these functions to convert 0, 0.5, 1, 1.5,
!   2.0, 2.5, ..., 8.0, 8.5, 9.0, 9.5 and 10.0 inch (resp., cm) to
!   cm (resp., inch).
!   ---------------------------------------------------------------

PROGRAM Conversion
   IMPLICIT NONE

     REAL, PARAMETER :: Initial = 0.0, Final = 10.0, Step = 0.5
     REAL            :: x

     x = Initial
     DO                        ! x = 0, 0.5, 1.0, ..., 9.0, 9.5, 10
        IF (x > Final) EXIT
        WRITE(*,*) x, 'cm = ',   Cm_to_Inch(x), 'inch and ', &
                         x, 'inch = ', Inch_to_Cm(x), 'cm'
       x = x + Step
    END DO

CONTAINS

! ---------------------------------------------------------------
! REAL FUNCTION Cm_to_Inch()
!    This function converts its real input in cm to inch.
! ---------------------------------------------------------------

    REAL FUNCTION Cm_to_Inch(cm)
       IMPLICIT NONE

        REAL, INTENT(IN) :: cm
        REAL, PARAMETER :: To_Inch = 0.3937           ! conversion factor

       Cm_to_Inch = To_Inch * cm
    END FUNCTION Cm_to_Inch

! ---------------------------------------------------------------
! REAL FUNCTION Inch_to_Cm()
!    This function converts its real input in inch to cm.
! ---------------------------------------------------------------

    REAL FUNCTION Inch_to_Cm(inch)
       IMPLICIT NONE

        REAL, INTENT(IN) :: inch
        REAL, PARAMETER :: To_Cm = 2.54               ! conversion factor

       Inch_to_Cm = To_Cm * inch
    END FUNCTION Inch_to_Cm

END PROGRAM Conversion
Click here to download this program.




PROGRAM INPUT AND OUTPUT
The following is the output from the above program.

 0.E+0cm = 0.E+0inch and 0.E+0inch = 0.E+0cm
 0.5cm = 0.196850002inch and 0.5inch = 1.26999998cm
 1.cm = 0.393700004inch and 1.inch = 2.53999996cm
 1.5cm = 0.590550005inch and 1.5inch = 3.80999994cm
 2.cm = 0.787400007inch and 2.inch = 5.07999992cm
 2.5cm = 0.984250009inch and 2.5inch = 6.3499999cm
 3.cm = 1.18110001inch and 3.inch = 7.61999989cm
 3.5cm = 1.37794995inch and 3.5inch = 8.88999939cm
 4.cm = 1.57480001inch and 4.inch = 10.1599998cm
 4.5cm = 1.77165008inch and 4.5inch = 11.4300003cm
 5.cm = 1.96850002inch and 5.inch = 12.6999998cm
 5.5cm = 2.16534996inch and 5.5inch = 13.9699993cm
 6.cm = 2.36220002inch and 6.inch = 15.2399998cm
 6.5cm = 2.55905008inch and 6.5inch = 16.5100002cm
 7.cm = 2.75589991inch and 7.inch =                  17.7799988cm
 7.5cm = 2.95274997inch and 7.5inch                  = 19.0499992cm
 8.cm = 3.14960003inch and 8.inch =                  20.3199997cm
 8.5cm = 3.34645009inch and 8.5inch                  = 21.5900002cm
 9.cm = 3.54330015inch and 9.inch =                  22.8600006cm
 9.5cm = 3.74014997inch and 9.5inch                  = 24.1299992cm
 10.cm = 3.93700004inch and 10.inch                  = 25.3999996cm


DISCUSSION

       Function Cm_to_Inch() converts its formal argument cm to inch, which is returned as the function value.
        Please note that the constant 0.3937 is defined as a PARAMETER.
       Function Inch_to_Cm() converts its formal argument inch to cm, which is returned as the function value.
        Please note that the constant 2.54 is defined as a PARAMETER.
       The main program uses DO-EXIT-END DO to generate 0, 0.5, 1, 1.5, ..., 8, 8.5, 9, 9.5 and 10. For each
        value, Cm_to_Inch() and Inch_to_Cm() are called to perform the desired conversion.


HERON'S FORMULA FOR COMPUTING TRIANGLE AREA - REVISITED

PROBLEM STATEMENT
Given a triangle with side lengths a, b and c, its area can be computed using the Heron's formula:




where s is the half of the perimeter length:




In order for a, b and c to form a triangle, two conditions must be satisfied. First, all side lengths
must be positive:




Second, the sum of any two side lengths must be greater than the third side length:




Write a program to read in three real values and use a function for testing the conditions and
another function for computing the area. Should the conditions fail, your program must keep
asking the user to re-enter the input until the input form a triangle. Then, the other function is
used to compute the area.
SOLUTION
! --------------------------------------------------------------------
!    This program uses Heron's formula to compute the area of a
! triangle. It "contains" the following functions;
!    (1) LOGICAL function TriangleTest() -
!         this function has three real formal arguments and tests
!         to see if they can form a triangle. If they do form a
!         triangle, this function returns .TRUE.; otherwise, it
!         returns .FALSE.
!    (2) REAL function TriangleArea() -
!         this functions has three real formal arguments considered
!         as three sides of a triangle and returns the area of this
!         triangle.
! --------------------------------------------------------------------

PROGRAM HeronFormula
   IMPLICIT NONE

     REAL :: a, b, c, TriangleArea

     DO
        WRITE(*,*) 'Three sides of a triangle please --> '
        READ(*,*)   a, b, c
        WRITE(*,*) 'Input sides are ', a, b, c
        IF (TriangleTest(a, b, c)) EXIT ! exit if not a triangle
        WRITE(*,*) 'Your input CANNOT form a triangle. Try again'
     END DO

     TriangleArea = Area(a, b, c)
     WRITE(*,*) 'Triangle area is ', TriangleArea

CONTAINS

!   --------------------------------------------------------------------
!   LOGICAL FUNCTION TriangleTest() :
!      This function receives three REAL numbers and tests if they form
!   a triangle by testing:
!      (1) all arguments must be positive, and
!      (2) the sum of any two is greater than the third
!   If the arguments form a triangle, this function returns .TRUE.;
!   otherwise, it returns .FALSE.
!   --------------------------------------------------------------------

     LOGICAL FUNCTION TriangleTest(a, b, c)
        IMPLICIT NONE

          REAL, INTENT(IN) :: a, b, c
          LOGICAL          :: test1, test2

        test1 = (a > 0.0) .AND. (b > 0.0) .AND. (c > 0.0)
        test2 = (a + b > c) .AND. (a + c > b) .AND. (b + c > a)
        TriangleTest = test1 .AND. test2 ! both must be .TRUE.
     END FUNCTION TriangleTest

! --------------------------------------------------------------------
! REAL FUNCTION Area() :
!    This function takes three real number that form a triangle, and
! computes and returns the area of this triangle using Heron's formula.
! --------------------------------------------------------------------

    REAL FUNCTION Area(a, b, c)
       IMPLICIT NONE

        REAL, INTENT(IN) :: a, b, c
        REAL             :: s

       s    = (a + b + c) / 2.0
       Area = SQRT(s*(s-a)*(s-b)*(s-c))
    END FUNCTION Area

END PROGRAM       HeronFormula

Click here to download this program.




PROGRAM INPUT AND OUTPUT
The following is the output from the above program.

Three sides of a triangle please -->
-3.0 4.0 5.0
Input sides are -3., 4., 5.
Your input CANNOT form a triangle. Try again
Three sides of a triangle please -->
1.0 3.0 4.0
Input sides are 1., 3., 4.
Your input CANNOT form a triangle. Try again
Three sides of a triangle please -->
6.0 8.0 10.0
Input sides are 6., 8., 10.
Triangle area is 24.


DISCUSSION

       LOGICAL function TriangleTest() receives three REAL values. The result of the first test condition is saved
        to a local LOGICAL variable test1, while the result of the second condition is saved to another LOGICAL
        variable test2. Since both conditions must be true to have a triangle, test1 and test2 are .AND.ed and the
        result goes into the function name so that it could be returned.
       REAL function Area is simple and does not require further discussion. However, please note that Area()
        has three formal arguments whose names are identical to the three global variables declared in the main
        program. By Scope Rule 3, they are different entities and do not cause any conflicts.
       The main program has a DO-EXIT-END DO loop. In each iteration, it asks for three real values. These
        values are sent to LOGICAL function TriangleTest() for testing. If the returned value is .TRUE., the input
        form a triangle and the control of execution exits. Then, the area is computed with REAL function Area().
        If the returned value is .FALSE., this function displays a message and goes back asking for a new set of
        values.


COMMUTING THE COMBINATORIAL COEFFICIENT
PROBLEM STATEMENT

The combinatorial coefficient C(n,r) is defined as follows:




where 0 <= r <= n must hold. Write a program that keeps reading in values for n and r, exits if
both values are zeros, uses a LOGICAL function to test if 0 <= r <= n holds, and computes
C(n,r) with an INTEGER function.


SOLUTION
!   ---------------------------------------------------------------
!   This program computes the combinatorial coefficient C(n,r):
!
!                         n!
!           C(n,r) = -------------
!                     r! x (n-r)!
!
!   It asks for two integers and uses Cnr(n,r) to compute the value.
!   If 0 <= r <= n does not hold, Cnr() returns -1 so that the main
!   program would know the input values are incorrect. Otherwise,
!   Cnr() returns the desired combinatorial coefficient.
!
!   Note that if the input values are zeros, this program stops.
!   ---------------------------------------------------------------

PROGRAM Combinatorial
   IMPLICIT  NONE

     INTEGER :: n, r, Answer

     DO
        WRITE(*,*)
        WRITE(*,*) "Two integers n and r (0 <= r <= n) please "
        WRITE(*,*) "0 0 to stop --> "
        READ(*,*)    n, r
        IF (n == 0 .AND. r == 0) EXIT
        WRITE(*,*) "Your input:"
        WRITE(*,*) " n        = ", n
        WRITE(*,*) " r        = ", r
        Answer = Cnr(n, r)
        IF (Answer < 0) THEN
            WRITE(*,*) "Incorrect input"
        ELSE
            WRITE(*,*) " C(n,r) = ", Answer
        END IF
     END DO
CONTAINS

!   ---------------------------------------------------------------
!   INTEGER FUNCTION Cnr(n,r)
!      This function receives n and r, uses LOGICAL function Test()
!   to verify if the condition 0 <= r <= n holds, and uses
!   Factorial() to compute n!, r! and (n-r)!.
!   ---------------------------------------------------------------

     INTEGER FUNCTION Cnr(n, r)
        IMPLICIT NONE
        INTEGER, INTENT(IN) :: n, r

        IF (Test(n,r)) THEN
           Cnr = Factorial(n)/(Factorial(r)*Factorial(n-r))
        ELSE
           Cnr = -1
        END IF
     END FUNCTION Cnr

!   ---------------------------------------------------------------
!   LOGICAL FUNCTION Test()
!      This function receives n and r. If 0 <= r <= n holds, it
!   returns .TRUE.; otherwise, it returns .FALSE.
!   ---------------------------------------------------------------

     LOGICAL FUNCTION Test(n, r)
        IMPLICIT NONE
        INTEGER, INTENT(IN) :: n, r

        Test = (0 <= r) .AND. (r <= n)
     END FUNCTION Test

!   ---------------------------------------------------------------
!   INTEGER FUNCTION Factorial()
!      This function receives a non-negative integer and computes
!   its factorial.
!   ---------------------------------------------------------------

     INTEGER FUNCTION Factorial(k)
        IMPLICIT NONE
        INTEGER, INTENT(IN) :: k
        INTEGER             :: Ans, i

        Ans = 1
        DO i = 1, k
           Ans = Ans * i
        END DO
        Factorial = Ans
     END FUNCTION Factorial

END PROGRAM Combinatorial
Click here to download this program.
PROGRAM INPUT AND OUTPUT
The following is the output from the above program.

Two integers n and r (0 <= r <= n) please
0 0 to stop -->
10 4
Your input:
  n      = 10
  r      = 4
  C(n,r) = 210

Two integers n and r (0 <= r <= n) please
0 0 to stop -->
7 6
Your input:
  n      = 7
  r      = 6
  C(n,r) = 7

Two integers n and r (0 <= r <= n) please
0 0 to stop -->
4 8
Your input:
  n      = 4
  r      = 8
Incorrect input

Two integers n and r (0 <= r <= n) please
0 0 to stop -->
-3 5
Your input:
  n      = -3
  r      = 5
Incorrect input

Two integers n and r (0 <= r <= n) please
0 0 to stop -->
0 0

In the sample output above, please note the error messages indicating that condition 0 <= r <= n does not hold.
Please also note that the program stops when the input values are 0 and 0.




DISCUSSION

       The function that computes the combinatorial coefficient is INTEGER function Cnr(n,r). Since it must
        compute n!, r! and (n-r)!, it would be better to write a function to compute the factorial of an integer. In
        this way, Cnr() would just use Factorial() three times rather than using three DO-loops. Note that if the
        condition does not hold, Cnr() returns -1.
       INTEGER function Factorial() computes the factorial of its formal argument.
       LOGICAL function Test() is very simple. If condition 0 <= r <= n holds, Test receives .TRUE.; otherwise, it
        receives .FALSE.
        The main program has a DO-EXIT-END DO. It asks for two integers and indicates that the program will
         stop if both values are zeros. Then, Cnr() is used for computing the combinatorial coefficient. If the
         returned value is negative, the input values do not meet the condition and an error message is displayed.
         Otherwise, the combinatorial coefficient is shown.
        Note that all three functions are internal functions of the main program.


COMPUTING SQUARE ROOTS WITH NEWTON'S METHOD

PROBLEM STATEMENT
We have discussed Newton's Method for computing the square root of a positive number. Let the given number
be b and let x be a rough guess of the square root of b. Newton's method suggests that a better guess, New x can
be computed as follows:




One can start with b as a rough guess and compute New x; from New x, one can generate a even
better guess, until two successive guesses are very close. Either one could be considered as the
square root of b.

Write a function MySqrt() that accepts a formal argument and uses Newton's method to
computes its square root. Then, write a main program that reads in an initial value, a final value,
and a step size, and computes the square roots of these successive values with Newton'e method
and Fortran's SQRT() function, and determines the absolute error.


SOLUTION
!   ---------------------------------------------------------------
!   This program contains a function MySqrt() that uses Newton's
!   method to find the square root of a positive number. This is
!   an iterative method and the program keeps generating better
!   approximation of the square root until two successive
!   approximations have a distance less than the specified tolerance.
!   ---------------------------------------------------------------

PROGRAM SquareRoot
   IMPLICIT NONE

     REAL       :: Begin, End, Step
     REAL       :: x, SQRTx, MySQRTx, Error

     READ(*,*)      Begin, End, Step                       ! read in init, final and step
     x = Begin                                             ! x starts with the init value
     DO
        IF (x >     End) EXIT                              !   exit if x >     the final value
        SQRTx       = SQRT(x)                              !   find square     root with SQRT()
        MySQRTx     = MySqrt(x)                            !   do the same     with my sqrt()
        Error       = ABS(SQRTx - MySQRTx)                 !   compute the     absolute error
        WRITE(*,*) x, SQRTx, MySQRTx, Error   ! display the results
        x = x + Step                     ! move on to the next value
     END DO

CONTAINS

!   ---------------------------------------------------------------
!   REAL FUNCTION MySqrt()
!      This function uses Newton's method to compute an approximate
!   of a positive number. If the input value is zero, then zero is
!   returned immediately. For convenience, the absolute value of
!   the input is used rather than kill the program when the input
!   is negative.
!   ---------------------------------------------------------------

     REAL FUNCTION MySqrt(Input)
        IMPLICIT NONE
        REAL, INTENT(IN) :: Input
        REAL             :: X, NewX
        REAL, PARAMETER :: Tolerance = 0.00001

        IF (Input == 0.0) THEN             ! if the input is zero
           MySqrt = 0.0                    !    returns zero
        ELSE                               ! otherwise,
           X = ABS(Input)                  !    use absolute value
           DO                              !    for each iteration
               NewX = 0.5*(X + Input/X)    !       compute a new approximation
               IF (ABS(X - NewX) < Tolerance) EXIT ! if very close, exit
               X = NewX                    !       otherwise, keep the new one
           END DO
           MySqrt = NewX
        END IF
     END FUNCTION MySqrt

END PROGRAM       SquareRoot
Click              here                to               download                this               program.



PROGRAM INPUT AND OUTPUT
If the input values of Begin, End and Step are 0.0, 10.0 and 0.5, the following is the output from the above
program.

0.E+0, 0.E+0, 0.E+0, 0.E+0
0.5, 0.707106769, 0.707106769, 0.E+0
1., 1., 1., 0.E+0
1.5, 1.22474492, 1.2247448, 1.192092896E-7
2., 1.41421354, 1.41421354, 0.E+0
2.5, 1.58113885, 1.58113885, 0.E+0
3., 1.73205078, 1.7320509, 1.192092896E-7
3.5, 1.87082875, 1.87082863, 1.192092896E-7
4., 2., 2., 0.E+0
4.5, 2.12132025, 2.12132025, 0.E+0
5., 2.23606801, 2.23606801, 0.E+0
5.5, 2.34520793, 2.34520793, 0.E+0
6., 2.44948983, 2.44948959, 2.384185791E-7
6.5, 2.54950976, 2.54950976, 0.E+0
7., 2.64575124, 2.64575148, 2.384185791E-7
7.5, 2.73861289, 2.73861265, 2.384185791E-7
8., 2.82842708, 2.82842708, 0.E+0
8.5, 2.91547585, 2.91547585, 0.E+0
9., 3., 3., 0.E+0
9.5, 3.08220696, 3.08220696, 0.E+0
10., 3.1622777, 3.1622777, 0.E+0


DISCUSSION
This program has nothing special. Please refer to the discussion of Newton's method for the computation details of
function MySqrt(). However, there is one thing worth to be mentioned. Since the formal argument Input is
declared with INTENT(IN), it cannot be changed in function MySqrt(). Therefore, the value of the formal argument
Input is copied to X and used in square root computation.


DESIGNING A GREATEST COMMON DIVISOR FUNCTION

PROBLEM STATEMENT
We have seen Greatest Common Divisor computation. This problem uses the same idea; but the computation is
performed                          with                            a                            function.



SOLUTION
!   ---------------------------------------------------------
!   This program computes the GCD of two positive integers
!   using the Euclid method. Given a and b, a >= b, the
!   Euclid method goes as follows: (1) dividing a by b yields
!   a reminder c; (2) if c is zero, b is the GCD; (3) if c is
!   no zero, b becomes a and c becomes c and go back to
!   Step (1). This process will continue until c is zero.
!
!   Euclid's algorithm is implemented as an INTEGER function
!   GCD().
!   ---------------------------------------------------------

PROGRAM GreatestCommonDivisor
   IMPLICIT NONE

     INTEGER      :: a, b

     WRITE(*,*) 'Two positive integers please --> '
     READ(*,*) a, b
     WRITE(*,*) 'The GCD of is ', GCD(a, b)

CONTAINS

! ---------------------------------------------------------
! INTEGER FUNCTION GCD():
!    This function receives two INTEGER arguments and
! computes their GCD.
! ---------------------------------------------------------

    INTEGER FUNCTION GCD(x, y)
       IMPLICIT NONE

        INTEGER, INTENT(IN) :: x, y                        ! we need x and y here
        INTEGER             :: a, b, c

        a = x                               !   if x <= y, swap x and y
        b = y                               !   since x and y are declared with
        IF (a <= b) THEN                    !   INTENT(IN), they cannot be
           c = a                            !   involved in this swapping process.
           a = b                            !   So, a, b and c are used instead.
           b = c
        END IF

        DO
           c = MOD(a, b)
           IF (c == 0) EXIT
           a = b
           b = c
        END DO

       GCD = b
    END FUNCTION         GCD

END PROGRAM        GreatestCommonDivisor

Click here to download this program.




PROGRAM INPUT AND OUTPUT
The following is the output from the above program.

Two positive integers please -->
46332 71162
The GCD of is 26


DISCUSSION
Nothing special is here. As in the previous example, since x and y are declared with INTENT(IN), their values cannot
be modified and therefore their values are copies to a and b to be used in other computation.




FINDING ALL PRIME NUMBERS IN THE RANGE OF 2 AND N - REVISITED

PROBLEM STATEMENT
We have discussed a method for finding all prime numbers in the range of 2 and N previously. In fact, one can
design a function with an INTEGER formal argument and returns .TRUE. if the argument is a prime number. Then,
it    is   used   to    test   if   the   integers   in   the   range   of   2   and   N   are   integers.



SOLUTION
!    --------------------------------------------------------------------
!    This program finds all prime numbers in the range of 2 and an
!    input integer.
!    --------------------------------------------------------------------

PROGRAM Primes
   IMPLICIT NONE

      INTEGER     :: Range, Number, Count

      Range = GetNumber()
      Count = 1                            ! input is correct. start counting
      WRITE(*,*)                           ! since 2 is a prime
      WRITE(*,*) 'Prime number #', Count, ': ', 2
      DO Number = 3, Range, 2              ! try all odd numbers 3, 5, 7, ...
         IF (Prime(Number)) THEN
             Count = Count + 1             ! yes, this Number is a prime
             WRITE(*,*) 'Prime number #', Count, ': ', Number
         END IF
      END DO

      WRITE(*,*)
      WRITE(*,*)       'There are ', Count, ' primes in the range of 2 and ', Range

CONTAINS

!    --------------------------------------------------------------------
!    INTEGER FUNCTION GetNumber()
!       This function does not require any formal argument. It keeps
!    asking the reader for an integer until the input value is greater
!    than or equal to 2.
!    --------------------------------------------------------------------

      INTEGER FUNCTION GetNumber()
         IMPLICIT NONE

           INTEGER :: Input

         WRITE(*,*) 'What is the range ? '
         DO                                ! keep trying to read a good input
            READ(*,*) Input                ! ask for an input integer
            IF (Input >= 2) EXIT           ! if it is GOOD, exit
            WRITE(*,*) 'The range value must be >= 2. Your input = ', Input
            WRITE(*,*) 'Please try again:'      ! otherwise, bug the user
         END DO
         GetNumber = Input
      END FUNCTION GetNumber

! --------------------------------------------------------------------
! LOGICAL FUNCTION Prime()
!    This function receives an INTEGER formal argument Number. If it
! is a prime number, .TRUE. is returned; otherwise, this function
! returns .FALSE.
! --------------------------------------------------------------------

    LOGICAL FUNCTION Prime(Number)
       IMPLICIT NONE

        INTEGER, INTENT(IN) :: Number
        INTEGER             :: Divisor

       IF (Number < 2) THEN
          Prime = .FALSE.
       ELSE IF (Number == 2) THEN
          Prime = .TRUE.
       ELSE IF (MOD(Number,2) == 0) THEN
          Prime = .FALSE.
       ELSE
          Divisor = 3
          DO
              IF (Divisor*Divisor>Number .OR. MOD(Number,Divisor)==0)   EXIT
              Divisor = Divisor + 2
          END DO
          Prime = Divisor*Divisor > Number
       END IF
    END FUNCTION Prime

END PROGRAM       Primes
Click              here                 to            download   this     program.



PROGRAM INPUT AND OUTPUT
The following is the output from the above program.

What is the range ?
-10
The range value must be >= 2.                Your input = -10
Please try again:
0
The range value must be >= 2.                Your input = 0
Please try again:
60

Prime   number    #1: 2
Prime   number    #2: 3
Prime   number    #3: 5
Prime   number    #4: 7
Prime   number    #5: 11
Prime   number    #6: 13
Prime   number    #7: 17
Prime   number    #8: 19
Prime   number    #9: 23
Prime   number    #10: 29
Prime   number    #11: 31
Prime   number    #12: 37
Prime   number    #13: 41
Prime   number    #14: 43
Prime   number    #15: 47
Prime number #16: 53
Prime number #17: 59

There are 17 primes in the range of 2 and 60


DISCUSSION

        Function GetNumber() has no formal arguments. It keeps asking the user to input an integer that is
         greater than or equal to two. The valid input is returned as the function value.
        The core part of this program is LOGICAL function Prime(). It receives an INTEGER formal argument
         Number and returns .TRUE. if Number is a prime number. For the working detail and logic of this
         function, please click here to bring you to the example discussed previously.


THE BISECTION (BOZANO) EQUATION SOLVER

PROBLEM STATEMENT

Given a continuous equation f(x)=0 and two values a and b (a < b), if f(a)*f(b) < 0 (i.e., f(a) and
f(b) have opposite signs), it can be proved that there exists a root of f(x)=0 between a and b.
More precisely, there exists a c, a <= c <= b, such that f(c)=0 holds.

This result provides us with a method for solving equations. If we take the midpoint of a and b,
c=(a+b)/2, and computes its function value f(c), we have the following cases:

    1.   If f(c) is very small (i.e., smaller than a tolerance value), then c can be considered as a root of f(x)=0 and
         we are done.
    2.   Otherwise, since f(a) and f(b) have opposite signs, the sign of f(c) is either identical to that of f(a) or that
         of f(b).
               o If the sign of f(c) is different from that of f(a), then since f(a) and f(c) have opposite signs, f(x)=0
                     has a root in the range of a and c.
               o If the sign of f(c) is different from that of f(b), then since f(b) and f(c) have opposite signs, f(x)=0
                     has a root in the range of c and b.
    3.   Therefore, if the original interval [a,b] is replaced with [a,c] in the first case or replaced with [c,b] in the
         second case, the length of the interval is reduced by half. If continue this process, after a number of steps,
         the length of the interval could be very small. However, since this small interval still contains a root of
         f(x)=0, we actually find an approximation of that root!

Write a program that contains two functions: (1) Funct() - the function f(x) and (2) Solve() - the equation solver
based on the above theory. Then, reads in a and b and uses function Solve() to find a root in the range of a and b.
Note that before calling Solve, your program should check if f(a)*f(b)<0 holds.

Your function Funct() is given below:
This function has a root near -0.89.

Since this process keeps dividing the intervals into two equal halves, it is usually referred to as
the bisection method. It is also known as Bozano's method.


SOLUTION
!   --------------------------------------------------------------------
!      This program solves equations with the Bisection Method. Given
!   a function f(x) = 0. The bisection method starts with two values,
!   a and b such that f(a) and f(b) have opposite signs. That is,
!   f(a)*f(b) < 0. Then, it is guaranteed that f(x)=0 has a root in
!   the range of a and b. This program reads in a and b (Left and Right
!   in this program) and find the root in [a,b].
!      In the following, function f() is REAL FUNCTION Funct() and
!   solve() is the function for solving the equation.
!   --------------------------------------------------------------------

PROGRAM Bisection
   IMPLICIT NONE

     REAL, PARAMETER :: Tolerance = 0.00001
     REAL            :: Left, fLeft
     REAL            :: Right, fRight
     REAL            :: Root

     WRITE(*,*)    'This program can solves equation F(x) = 0'
     WRITE(*,*)    'Please enter two values Left and Right such that '
     WRITE(*,*)    'F(Left) and F(Right) have opposite signs.'
     WRITE(*,*)
     WRITE(*,*)    'Left and Right please --> '
     READ(*,*)     Left, Right              ! read in Left and Right

     fLeft = Funct(Left)                  ! compute their function values
     fRight = Funct(Right)
     WRITE(*,*)
     WRITE(*,*) 'Left = ', Left, '     f(Left) = ', fLeft
     WRITE(*,*) 'Right = ', Right, '    f(Right) = ', fRight
     WRITE(*,*)
     IF (fLeft*fRight > 0.0) THEN
        WRITE(*,*) '*** ERROR: f(Left)*f(Right) must be negative ***'
     ELSE
        Root = Solve(Left, Right, Tolerance)
        WRITE(*,*) 'A root is ', Root
     END IF

CONTAINS

!   --------------------------------------------------------------------
!   REAL FUNCTION Funct()
!      This is for function f(x). It takes a REAL formal argument and
!   returns the value of f() at x. The following is sample function
!   with a root in the range of -10.0 and 0.0. You can change the
!   expression with your own function.
!   --------------------------------------------------------------------
    REAL FUNCTION Funct(x)
       IMPLICIT NONE
       REAL, INTENT(IN) :: x
       REAL, PARAMETER :: PI = 3.1415926
       REAL, PARAMETER :: a = 0.8475

       Funct = SQRT(PI/2.0)*EXP(a*x) + x/(a*a + x*x)

    END FUNCTION   Funct

!   --------------------------------------------------------------------
!   REAL FUNCTION Solve()
!      This function takes Left - the left end, Right - the right end,
!   and Tolerance - a tolerance value such that f(Left)*f(Right) < 0
!   and find a root in the range of Left and Right.
!      This function works as follows. Because of INTENT(IN), this
!   function cannot change the values of Left and Right and therefore
!   the values of Left and Right are saved to a and b.
!      Then, the middle point c=(a+b)/2 and its function value f(c)
!   is computed. If f(a)*f(c) < 0, then a root is in [a,c]; otherwise,
!   a root is in [c,b]. In the former case, replacing b and f(b) with
!   c and f(c), we still maintain that a root in [a,b]. In the latter,
!   replacing a and f(a) with c and f(c) will keep a root in [a,b].
!   This process will continue until |f(c)| is less than Tolerance and
!   hence c can be considered as a root.
!   --------------------------------------------------------------------

    REAL FUNCTION Solve(Left, Right, Tolerance)
       IMPLICIT NONE
       REAL, INTENT(IN) :: Left, Right, Tolerance
       REAL             :: a, Fa, b, Fb, c, Fc

       a = Left                            ! save Left and Right
       b = Right

       Fa = Funct(a)                     !   compute the function values
       Fb = Funct(b)
       IF (ABS(Fa) < Tolerance) THEN     !   if f(a) is already small
          Solve = a                      !   then a is a root
       ELSE IF (ABS(Fb) < Tolerance) THEN       ! is f(b) is small
          Solve = b                      !   then b is a root
       ELSE                              !   otherwise,
          DO                             !   iterate ....
             c = (a + b)/2.0             !     compute the middle point
             Fc = Funct(c)               !     and its function value
             IF (ABS(Fc) < Tolerance) THEN      ! is it very small?
                 Solve = c               !   yes, c is a root
                 EXIT
             ELSE IF (Fa*Fc < 0.0) THEN !    do f(a)*f(c) < 0 ?
                 b = c                   !   replace b with c
                 Fb = Fc                 !   and f(b) with f(c)
             ELSE                        !   then f(c)*f(b) < 0 holds
                 a = c                   !   replace a with c
                 Fa = Fc                 !   and f(a) with f(c)
             END IF
          END DO                         !   go back and do it again
       END IF
    END FUNCTION        Solve

END PROGRAM       Bisection
Click              here                  to                download                   this                program.



PROGRAM INPUT AND OUTPUT
The following is the output from the above program with correct input:

This program solves equation F(x) = 0
Please enter two values Left and Right such that
F(Left) and F(Right) have opposite signs.

Left and Right please -->
-10.0 0.0

Left = -10.          f(Left) = -9.902540594E-2
Right = 0.E+0         f(Right) = 1.25331414

A root is -0.89050293
The following output shows that the function values of the input do not have opposite signs and hence program
stops.

This program solves equation F(x) = 0
Please enter two values Left and Right such that
F(Left) and F(Right) have opposite signs.

Left and Right please -->
-10.0 -1.0

Left = -10.          f(Left) = -9.902540594E-2
Right = -1.         f(Right) = -4.495930672E-2

*** ERROR: f(Left)*f(Right) must be negative ***


DISCUSSION

       Function Funct(x) returns the function value at x.
       The heart of this program is function Solve(). It takes three arguments Left, Right and Tolerance and finds
        a root in the range of Left and Right.
       Since arguments Left and Right are declared with INTENT(IN), their values cannot be changed. As a result,
        their values are copied to variables a and b.
       The function values of a and b are stored in Fa and Fb.
       If the absolute value of Fa (resp., Fb) is very small, one can consider a (resp., b) as a root of the given
        equation.
       Otherwise, the midpoint c of a and b and its function value Fc are computed. If Fc is very small, c is
        considered as a root.
       Then if Fa and Fc have opposite signs, replacing b and Fb with c and Fc, respectively.
       If Fa and Fc have the same sign, then Fc and Fb must have the opposite sign. In this case, a and Fa are
        replaced with c and Fc.
        Either way, the original interval [a,b] is replaced with a new one with half length. This process continues
         until the absolute value of the function value at c, Fc, is smaller than the given tolerance value. In this
         case, c is considered a root of the given equation.



A Brief Introduction to Modules
WHAT IS A MODULE?


In many previous example, there are several internal functions packed at the end of the main
program. In fact, many of them such as Factorial(), Combinatorial(), GCD(), and Prime() are
general functions that can also be used in other programs. To provide the programmers with a
way of packing commonly used functions into a few tool-boxes, Fortran 90 has a new capability
called modules. It has a syntactic form very similar to a main program, except for something that
are specific to modules.

SYNTAX
The following is the syntax of a module:

MODULE module-name
   IMPLICIT NONE
   [specification part]
CONTAINS
   [internal-functions]
END MODULE module-name

The differences between a program and a module are the following:

        The structure of a module is almost identical to the structure of a program.
        However, a module starts with the keyword MODULE and ends with END MODULE rather than
         PROGRAM and END PROGRAM.
        A module has specification part and could contain internal function; but, it does not have any statements
         between the specification part and the keyword CONTAINS.
        Consequently, a module does not contains statements to be executed as in a program. A module can only
         contain declarations and functions to be used by other modules and programs. This is perhaps one of the
         most important difference. As a result, a module cannot exist alone; it must be used with other modules
         and a main program



SHORT EXAMPLES
Here are some short examples of modules:

        The following is a very simple module. It has the specification part and does not have any internal
         function. The specification part has two REAL PARAMETERs, namely PI and g and one INTEGER variable
         Counter.
      MODULE SomeConstants
         IMPLICIT NONE
         REAL, PARAMETER :: PI = 3.1415926
         REAL, PARAMETER :: g = 980
         INTEGER         :: Counter
      END MODULE SomeConstants

     The following module SumAverage does not have any specification part (hence no IMPLICIT NONE); but it
      does contain two functions, Sum() and Average(). Please note that function Average() uses Sum() to
      compute the sum of three REAL numbers.

      MODULE     SumAverage

      CONTAINS
         REAL FUNCTION Sum(a, b, c)
            IMPLICIT NONE
            REAL, INTENT(IN) :: a, b, c
            Sum = a + b + c
         END FUNCTION Sum

          REAL FUNCTION Average(a, b, c)
             IMPLICIT NONE
             REAL, INTENT(IN) :: a, b, c
             Average = Sum(a,b,c)/2.0
          END FUNCTION Average

      END MODULE       SumAverage

     The following module DegreeRadianConversion contains two PARAMETERs, PI and Degree180, and two
      functions DegreeToRadian() and RadianToDegree(). The two PARAMETERs are used in the conversion
      functions. Note that these parameters are global to the two functions. See scope rules for the details.

      MODULE DegreeRadianConversion
         IMPLICIT NONE
         REAL, PARAMETER :: PI = 3.1415926
         REAL, PARAMETER :: Degree180 = 180.0

          REAL FUNCTION DegreeToRadian(Degree)
             IMPLICIT NONE
             REAL, INTENT(IN) :: Degree
             DegreeToRadian = Degree*PI/Degree180
          END FUNCTION DegreeToRadian

          REAL FUNCTION RadianToDegree(radian)
             IMPLICIT NONE
             REAL, INTENT(IN) :: Radian
             RadianToDegree = Radian*Degree180/PI
          END FUNCTION RadianToDegree

      END MODULE       DegreeRadianConversion



HOW TO USE A MODULE?
Once a module is written, its global entities (i.e., global PARAMETERs, global variables and
internal functions) can be made available to other modules and programs. The program or
module that wants to use a particular module must have a USE statement at its very beginning.
The USE statement has one of the following forms:

SYNTAX
The following form is the syntax of a module:

USE     module-name

USE     module-name, ONLY: name-1, name-2, ..., name-n

Here is the meaning of these USEs:

      The first indicates that the current program or module wants to use the module whose name is module-name.
      For example, the following main program indicates that it wants to use the content of module
      SomeConstants:

      PROGRAM MainProgram
         USE   SomeConstants
         IMPLICIT NONE
         ..........
      END PROGRAM MainProgram

      Once a USE is specified in a program or in a module, every global entities of that used module (i.e.,
      PARAMETERs, variables, and internal functions) becomes available to this program or this module. For
      example, if module SomeConstants is:

      MODULE SomeConstants
         IMPLICIT NONE
         REAL, PARAMETER :: PI = 3.1415926
         REAL, PARAMETER :: g = 980
         INTEGER         :: Counter
      END MODULE SomeConstants

      Then, program MainProgram can access to PARAMETERs PI and g and variable Counter.

      However, under many circumstances, a program or a module does not want to use everything of the used
      module. For example, if program MainProgram only wants to use PARAMETER PI and variable Counter and
      does not want to use g, then the second form becomes very useful. In this case, one should add the keyword
      ONLY followed by a colon :, followed by a list of names that the current program or module wants to use.

      PROGRAM MainProgram
         USE   SomeConstants, ONLY: PI, Counter
         IMPLICIT NONE
         ..........
      END PROGRAM MainProgram

          Thus, MainProgram can use PI and Counter of module SomeConstants; but, MainProgram cannot use g!
    USE with ONLY: is very handy, because it could only "import" those important and
    vital information and "ignore" those un-wanted ones.

   There is a third, not recommended form, that can help to rename the names of a module locally. If a
    module has a name abc, in a program with a USE, one can use a new name for abc. The way of writing this
    renaming is the following:
     new-name => name-in-module

    For example, if one wants to give a new name to Counter of module SomeConstants, one should write:

    NewCounter => Counter

    Thus, in a program, whenever it uses MyCounter, this program is actually using Counter of module
    SomeConstants.

    This kind of renaming can be used with ONLY: or without ONLY:.

    PROGRAM MainProgram
       USE   SomeConstants, ONLY: PI, MyCounter => Counter
       IMPLICIT NONE
       ..........
    END PROGRAM MainProgram

    The above example wants to use PI of module SomeConstants. In this case, when program MainProgram
    uses PI, it actually uses the PI of module SomeConstants. Program MainProgram also uses Counter of
    module SomeConstants; but, in this case, since there is a renaming, when MainProgram uses MyCounter
    it actually uses Counter of module SomeConstants.

    Renaming does not require ONLY:. In the following, MainProgram can use all contents
    of module SomeConstants. When MainProgram uses PI and g, it uses the PI and g of
    module SomeConstants; however, when MainProgram uses MyCounter, it actually
    uses Counter of module SomeConstants.

    PROGRAM MainProgram
       USE   SomeConstants, MyCounter => Counter
       IMPLICIT NONE
       ..........
    END PROGRAM MainProgram

    Why do we need renaming? It is simple. In your program, you may have a variable whose name is
    identical to an entity of a module that is being USEed. In this case, renaming the variable name would
    clear the ambiguity.

    PROGRAM MainProgram
       USE   SomeConstants, GravityConstant => g
       IMPLICIT NONE
       INTEGER :: e, f, g
       ..........
    END PROGRAM MainProgram
        In the above example, since MainProgram has a variable called g, which is the same as PARAMETER g in
        module SomeConstants. By renaming g of SomeConstants, MainProgram can use variable g for the
        variable and GravityConstant for the PARAMETER g in module SomeConstants.

        However, renaming is not a recommended feature. You should avoid using it
        whenever possible.




COMPILE PROGRAMS WITH MODULES


Normally, your programs and modules are stored in different files, all with filename suffix .f90.
When you compile your program, all involved modules must also be compiled. For example, if
your program is stored in a file main.f90 and you expect to use a few modules stored in files
compute.f90, convert.f90 and constants.f90, then the following command will compile your
program and all three modules, and make an executable file a.out

f90 compute.f90 convert.f90 constants.f90 main.f90
If you do not want the executable to be called a.out, you can do the following:

f90 compute.f90 convert.f90 constants.f90 main.f90 -o main
In this case, -o main instructs the Fortran compiler to generate an executable main instead of a.out.

Different compilers may have different "personalities." This means the way of compiling
modules and your programs may be different from compilers to compilers. If you have
several modules, say A, B, C, D and E, and C uses A, D uses B, and E uses A, C and D, then
the safest way to compile your program is the following command:

f90 A.f90 B.f90 C.f90 D.f90 E.f90 main.f90
That is, list those modules that do not use any other modules first, followed by those modules that only use
those listed modules, followed by your main program.

You may also compile your modules separately. For example, you may want to write and
compile all modules before you start working on your main program. The command for you to
compile a single program or module without generating an executable is the following:

f90 -c compute.f90
where -c means "compile the program without generating an executable." The output from Fortran compiler is a
file whose name is identical to the program file's name but with a file extension .o. Therefore, the above command
generates a file compute.o. The following commands compile all of your modules:

f90 -c compute.f90
f90 -c convert.f90
f90 -c constants.f90
After compiling these modules, your will see at least three files, compute.o, convert.o and constants.o.
After completing your main program, you have two choices: (1) compile your main program
separately, or (2) compile your main with all "compiled" modules. In the former, you perhaps
use the following command to compile your main program:

f90 -c main.f90
Then, you will have a file main.o. Now you have four .o files main.o, compute.o, convert.o and constants.o. To
pull them into an executable, you need

f90 compute.o convert.o constants.o main.o
This would generate a.out. Or, you may want to use

f90 compute.o convert.o constants.o main.o -o main
so that your executable could be called main.

If you want compile your main program main.f90 with all .o files, then you need

f90 compute.o convert.o constants.o main.c
or

f90 compute.o convert.o constants.o main.c -o main

Note that the order of listing these .o files may be important. Please consult the rules
discussed earlier.




FACTORIAL AND COMBINATORIAL COEFFICIENT

PROBLEM STATEMENT

The combinatorial coefficient C(n,r) is defined as follows:




where 0 <= r <= n must hold. Write a module that contains two functions: (1) Factorial() and
(2) Combinatorial(). The former computes the factorial of its argument, while the latter uses the
former to compute the combinatorial coefficient. Then, write a main program that uses this
module.


SOLUTION
The following is the desired module:

! --------------------------------------------------------------------
! MODULE FactorialModule
!      This module contains two procedures: Factorial(n) and
!   Combinatorial(n,r). The first computes the factorial of an integer
!   n and the second computes the combinatorial coefficient of two
!   integers n and r.
!   --------------------------------------------------------------------

MODULE FactorialModule
   IMPLICIT NONE

CONTAINS

!   --------------------------------------------------------------------
!   FUNCTION Factorial() :
!      This function accepts a non-negative integers and returns its
!   Factorial.
!   --------------------------------------------------------------------

     INTEGER FUNCTION Factorial(n)
        IMPLICIT NONE

        INTEGER, INTENT(IN) :: n            ! the argument
        INTEGER             :: Fact, i      ! result

        Fact = 1                            ! initially, n!=1
        DO i = 1, n                         ! this loop multiplies
           Fact = Fact * i                  ! i to n!
        END DO
        Factorial = Fact

     END FUNCTION     Factorial

!   --------------------------------------------------------------------
!   FUNCTION Combinarotial():
!      This function computes the combinatorial coefficient C(n,r).
!   If 0 <= r <= n, this function returns C(n,r), which is computed as
!   C(n,r) = n!/(r!*(n-r)!). Otherwise, it returns 0, indicating an
!   error has occurred.
!   --------------------------------------------------------------------

     INTEGER FUNCTION Combinatorial(n, r)
        IMPLICIT NONE

        INTEGER, INTENT(IN) :: n, r
        INTEGER             :: Cnr

        IF (0 <= r .AND. r <= n) THEN     ! valid arguments ?
           Cnr = Factorial(n) / (Factorial(r)*Factorial(n-r))
        ELSE                              ! no,
           Cnr = 0                        ! zero is returned
        END IF
        Combinatorial = Cnr

     END FUNCTION     Combinatorial

END MODULE FactorialModule
Click here to download this program.
Here is the main program:

!   --------------------------------------------------------------------
!   PROGRAM ComputeFactorial:
!      This program uses MODULE FactorialModule for computing factorial
!   and combinatorial coefficients.
!   --------------------------------------------------------------------

PROGRAM      ComputeFactorial
   USE           FactorialModule                            ! use a module

     IMPLICIT      NONE

     INTEGER :: N, R

     WRITE(*,*)       'Two non-negative integers --> '
     READ(*,*)        N, R

     WRITE(*,*)       N,      '! = ', Factorial(N)
     WRITE(*,*)       R,      '! = ', Factorial(R)

     IF (R <= N) THEN                    ! if r <= n, do C(n,r)
        WRITE(*,*) 'C(', N, ',', R, ') = ', Combinatorial(N, R)
     ELSE                                ! otherwise, do C(r,n)
        WRITE(*,*) 'C(', R, ',', N, ') = ', Combinatorial(R, N)
     END IF

END PROGRAM        ComputeFactorial
Click               here          to                         download                  this                program.



PROGRAM INPUT AND OUTPUT
The following is the output from the above program.

Two non-negative integers -->
13 4
13! = 1932053504
4! = 24
C(13,4) = 221


DISCUSSION

        The computation of combinatorial coefficients has been discussed in an programming example, where
         functions Cnr(n,r) and Factorial(k) are internal functions of the main program.
        In this version, functions Factorial(n) and Combinatorial(n,r) are moved to a module called
         FactorialModule as internal functions of that module.
        Factorial(n) takes a non-negative integer and returns its factorial.
        Combinatorial(n,r) takes two non-negative integers n and r. If 0 <= r <= n, the combinatorial coefficient
         C(n,r) is returned; otherwise, 0 is returned.
        Note that in module FactorialModule, there is no variables global to its internal functions. All internal
         functions use their own internal (or local) variables.
       This module does not perform many checks as in a previous programming example. But, it is not difficult
        to add these tests.
       After moving the computation functions to a module, the main program becomes simpler. In the
        beginning, the main program must USES FactorialModule so that functions Factorial() and
        Combinatorial() can be accessed from within the main program.
       The main program reads in values for n and r. If r <= n, the combinatorial coefficient C(n,r) is computed by
        calling Combinatorial(n,r); otherwise, the main program computes Combinatorial(r,n).
       If the main program and module FactorialModule are stored in files fact-1p.f90 and fact-m.f90,
        respectively, then you can compile them together with the following command:
           f90 fact-m.f90 fact-1p.90

        or with the following that generates an executable called fact1p:

        f90 fact-m.f90 fact-1p.90 -o fact-1p


TRIGONOMETRIC FUNCTIONS USING DEGREE

PROBLEM STATEMENT

Trigonometric functions (i.e., sin(x) and cos(x)) use radian for their argument. Using modules,
you can design your own trigonometric functions that use degree. Write a module that contains
functions for converting radian to degree and degree to radian and sin(x) and cos(x) with
arguments in degree rather than in radian.


SOLUTION
The following is the desired module:

! --------------------------------------------------------------------
! MODULE MyTrigonometricFunctions:
!    This module provides the following functions and constants
!    (1) RadianToDegree()     - converts its argument in radian to
!                               degree
!    (2) DegreeToRadian()     - converts its argument in degree to
!                               radian
!    (3) MySIN()              - compute the sine of its argument in
!                               degree
!    (4) MyCOS()              - compute the cosine of its argument
!                               in degree
! --------------------------------------------------------------------

MODULE MyTrigonometricFunctions
   IMPLICIT  NONE

    REAL,    PARAMETER     ::   PI             =   3.1415926                ! some constants
    REAL,    PARAMETER     ::   Degree180      =   180.0
    REAL,    PARAMETER     ::   R_to_D         =   Degree180/PI
    REAL,    PARAMETER     ::   D_to_R         =   PI/Degree180

CONTAINS
!   --------------------------------------------------------------------
!   FUNCTION RadianToDegree():
!      This function takes a REAL argument in radian and converts it to
!   the equivalent degree.
!   --------------------------------------------------------------------

     REAL FUNCTION RadianToDegree(Radian)
        IMPLICIT NONE
        REAL, INTENT(IN) :: Radian

        RadianToDegree = Radian * R_to_D
     END FUNCTION RadianToDegree

!   --------------------------------------------------------------------
!   FUNCTION DegreeToRadian():
!      This function takes a REAL argument in degree and converts it to
!   the equivalent radian.
!   --------------------------------------------------------------------

     REAL FUNCTION DegreeToRadian(Degree)
        IMPLICIT NONE
        REAL, INTENT(IN) :: Degree

        DegreeToRadian = Degree * D_to_R
     END FUNCTION DegreeToRadian

!   --------------------------------------------------------------------
!   FUNCTION MySIN():
!      This function takes a REAL argument in degree and computes its
!   sine value. It does the computation by converting its argument to
!   radian and uses Fortran's sin().
!   --------------------------------------------------------------------

     REAL FUNCTION MySIN(x)
        IMPLICIT NONE
        REAL, INTENT(IN) :: x

        MySIN = SIN(DegreeToRadian(x))
     END FUNCTION MySIN

!   --------------------------------------------------------------------
!   FUNCTION MySIN():
!      This function takes a REAL argument in degree and computes its
!   cosine value. It does the computation by converting its argument to
!   radian and uses Fortran's cos().
!   --------------------------------------------------------------------

     REAL FUNCTION MyCOS(x)
        IMPLICIT NONE
        REAL, INTENT(IN) :: x

        MyCOS = COS(DegreeToRadian(x))
     END FUNCTION MyCOS

END MODULE MyTrigonometricFunctions
Click here to download this program.
Here is the main program:

!   -----------------------------------------------------------------------
!   PROGRAM TrigonFunctTest:
!      This program tests the functions in module MyTrigonometricFunctions.
!   Module MyTrigonometricFunctions is stored in file trigon.f90.
!   Functions in that module use degree rather than radian. This program
!   displays the sin(x) and cos(x) values for x=-180, -170, ..., 0, 10, 20,
!   30, ..., 160, 170 and 180. Note that the sin() and cos() function
!   in module MyTrigonometricFunctions are named MySIN(x) and MyCOS(x).
!   -----------------------------------------------------------------------

PROGRAM TrigonFunctTest
   USE MyTrigonometricFunctions                       ! use a module

        IMPLICIT      NONE

        REAL   ::   Begin = -180.0                    ! initial value
        REAL   ::   Final = 180.0                     ! final value
        REAL   ::   Step =    10.0                    ! step size
        REAL   ::   x

        WRITE(*,*) 'Value of PI = ', PI
        WRITE(*,*)
        x = Begin                           ! start with 180 degree
        DO
           IF (x > Final) EXIT              ! if x > 180 degree, EXIT
           WRITE(*,*) 'x = ', x, 'deg    sin(x) = ', MySIN(x), &
                        ' cos(x) = ', MyCOS(x)
           x = x + Step                     ! advance x
        END DO

END PROGRAM           TrigonFunctTest
Click                  here           to              download          this   program.



PROGRAM INPUT AND OUTPUT
The following is the output from the above program.

Value of PI = 3.1415925

x   =    -180.deg       sin(x) = 8.742277657E-8    cos(x) = -1.
x   =    -170.deg       sin(x) = -0.173648298    cos(x) = -0.98480773
x   =    -160.deg       sin(x) = -0.342020214    cos(x) = -0.939692616
x   =    -150.deg       sin(x) = -0.50000006    cos(x) = -0.866025388
x   =    -140.deg       sin(x) = -0.642787635    cos(x) = -0.766044438
x   =    -130.deg       sin(x) = -0.766044438    cos(x) = -0.642787635
x   =    -120.deg       sin(x) = -0.866025388    cos(x) = -0.50000006
x   =    -110.deg       sin(x) = -0.939692616    cos(x) = -0.342020124
x   =    -100.deg       sin(x) = -0.98480773    cos(x) = -0.173648193
x   =    -90.deg       sin(x) = -1.   cos(x) = -4.371138829E-8
x   =    -80.deg       sin(x) = -0.98480773   cos(x) = 0.173648223
x   =    -70.deg       sin(x) = -0.939692616    cos(x) = 0.342020154
x   =    -60.deg       sin(x) = -0.866025448    cos(x) = 0.49999997
x   =    -50.deg       sin(x) = -0.766044438    cos(x) = 0.642787635
x   =       -40.deg      sin(x) = -0.642787576    cos(x) = 0.766044438
x   =       -30.deg      sin(x) = -0.5   cos(x) = 0.866025388
x   =       -20.deg      sin(x) = -0.342020124    cos(x) = 0.939692616
x   =       -10.deg      sin(x) = -0.173648179    cos(x) = 0.98480773
x   =       0.E+0deg      sin(x) = 0.E+0    cos(x) = 1.
x   =       10.deg      sin(x) = 0.173648179    cos(x) = 0.98480773
x   =       20.deg      sin(x) = 0.342020124    cos(x) = 0.939692616
x   =       30.deg      sin(x) = 0.5   cos(x) = 0.866025388
x   =       40.deg      sin(x) = 0.642787576    cos(x) = 0.766044438
x   =       50.deg      sin(x) = 0.766044438    cos(x) = 0.642787635
x   =       60.deg      sin(x) = 0.866025448    cos(x) = 0.49999997
x   =       70.deg      sin(x) = 0.939692616    cos(x) = 0.342020154
x   =       80.deg      sin(x) = 0.98480773    cos(x) = 0.173648223
x   =       90.deg      sin(x) = 1.   cos(x) = -4.371138829E-8
x   =       100.deg      sin(x) = 0.98480773    cos(x) = -0.173648193
x   =       110.deg      sin(x) = 0.939692616    cos(x) = -0.342020124
x   =       120.deg      sin(x) = 0.866025388    cos(x) = -0.50000006
x   =       130.deg      sin(x) = 0.766044438    cos(x) = -0.642787635
x   =       140.deg      sin(x) = 0.642787635    cos(x) = -0.766044438
x   =       150.deg      sin(x) = 0.50000006    cos(x) = -0.866025388
x   =       160.deg      sin(x) = 0.342020214    cos(x) = -0.939692616
x   =       170.deg      sin(x) = 0.173648298    cos(x) = -0.98480773
x   =       180.deg      sin(x) = -8.742277657E-8    cos(x) = -1.


DISCUSSION

             Module MyTrigonometricFunctions defines four constants:
                   1. PI,
                   2. Degree180,
                   3. R_to_D for radian to degree conversion, and
                   4. D_to_R for degree to radian conversion.
             Module MyTrigonometricFunctions contains four functions:
                   1. RadianToDegree() converts its radian argument to degree;
                   2. DegreeToRadian() converts its degree argument to radian;
                   3. MySIN() takes a degree argument and returns its sine value; and
                   4. MyCOS() takes a degree argument and returns its cosine value.
             In the module, MySIN(x) first converts its degree argument x to radian using function DegreeToRadian()
              and supplies the result to the Fortran SIN() function for computing the value of sine.
             The module uses MySIN() and MyCOS() rather than the same Fortran names SIN() and COS() to avoid
              confusion.
             As described in the use of modules, the main program can use variables, PARAMETERs and functions
              declared and defined in a module. Thus, the main program can use the value of PARAMETER PI defined in
              module MyTrigonometricFunctions. This is why there is no declaration of PI in the main program.
             In the main program, REAL variable x starts with 180 (in degree) and steps through 180 (in degree) with a
              step size 10. For each value of x, its SIN() and COS() is computed.
             If the main program and module MyTrigonometricFunctions are stored in files trigon.f90 and tri-test.f90,
              respectively, then you can compile them together with the following command:
                 f90 trigon.f90 tri-test.90

              or with the following that generates an executable called tri-test:

              f90 trigon.f90 tri-test.90 -o tri-test
A LITTLE PRIVACY - PUBLIC/PRIVATE


All global entities of a module, by default, can be accessed by a program or another module
using the USE statement. But, it is possible to set some restrictions that some entities are private.
A private entity of a module can only be accessed within that module. On the other hand, one
can explicitly list those entities that can be accessed from outside. This is done with the PUBLIC
and PRIVATE statements:

SYNTAX
The following is the syntax of a module:

PUBLIC     :: name-1, name-2, ..., name-n

PRIVATE :: name-1, name-2, ..., name-n

All entities listed in PRIVATE will not be accessible from outside of the module and all entities
listed in PUBLIC can be accessed from outside of the module. All not listed entities, by default,
can be accessed from outside of the module.

You can have many PRIVATE and PUBLIC statements.

In the following code segment, since VolumeOfDeathStar, SecretConstant and BlackKnight
are listed in a statement, they can only be used with the module. On the other hand, SkyWalker
and Princess are listed in PUBLIC, they can be accessed from outside of the module. There are
entities not listed: function WeaponPower() and DeathStar. By default, they are public and can
be accessed from outside of the module.

MODULE TheForce
   IMPLICIT  NONE

    INTEGER :: SkyWalker, Princess
    REAL    :: BlackKnight
    LOGICAL :: DeathStar

    REAL, PARAMETER :: SecretConstant = 0.123456

    PUBLIC :: SkyWalker, Princess
    PRIVATE :: VolumeOfDeathStar
    PRIVATE :: SecretConstant, BlackKnight

CONTAINS
   INTEGER FUNCTION VolumeOfDeathStar()
      ..........
   END FUNCTION WolumeOfDeathStar

    REAL FUNCTION         WeaponPower(SomeWeapon)
       ..........
    END FUNCTION
      ..........
END MODULE TheForce


A PROGRAMMING EXAMPLE
In a previous example of using degree in trigonometric functions, four constants and four functions are defined.
But, most of them are used in and meaningful to the module MyTrigonometricFunctions. Thus, one can make
them private so that they cannot be accessed from outside of this module. Here is a rewritten version:

! --------------------------------------------------------------------
! MODULE MyTrigonometricFunctions:
!    This module provides the following functions and constants
!    (1) RadianToDegree()     - converts its argument in radian to
!                               degree
!    (2) DegreeToRadian()     - converts its argument in degree to
!                               radian
!    (3) MySIN()              - compute the sine of its argument in
!                               degree
!    (4) MyCOS()              - compute the cosine of its argument
!                               in degree
! --------------------------------------------------------------------

MODULE MyTrigonometricFunctions
   IMPLICIT  NONE

     REAL,   PARAMETER     ::   PI           =   3.1415926             ! some constants
     REAL,   PARAMETER     ::   Degree180    =   180.0
     REAL,   PARAMETER     ::   R_to_D       =   Degree180/PI
     REAL,   PARAMETER     ::   D_to_R       =   PI/Degree180

     PRIVATE               :: Degree180, R_to_D, D_to_R
     PRIVATE               :: RadianToDegree, DegreeToRadian
     PUBLIC                :: MySIN, MyCOS

CONTAINS

!   --------------------------------------------------------------------
!   FUNCTION RadianToDegree():
!      This function takes a REAL argument in radian and converts it to
!   the equivalent degree.
!   --------------------------------------------------------------------

     REAL FUNCTION RadianToDegree(Radian)
        IMPLICIT NONE
        REAL, INTENT(IN) :: Radian

        RadianToDegree = Radian * R_to_D
     END FUNCTION RadianToDegree

!   --------------------------------------------------------------------
!   FUNCTION DegreeToRadian():
!      This function takes a REAL argument in degree and converts it to
!   the equivalent radian.
!   --------------------------------------------------------------------
     REAL FUNCTION DegreeToRadian(Degree)
        IMPLICIT NONE
        REAL, INTENT(IN) :: Degree

        DegreeToRadian = Degree * D_to_R
     END FUNCTION DegreeToRadian

!   --------------------------------------------------------------------
!   FUNCTION MySIN():
!      This function takes a REAL argument in degree and computes its
!   sine value. It does the computation by converting its argument to
!   radian and uses Fortran's sin().
!   --------------------------------------------------------------------

     REAL FUNCTION MySIN(x)
        IMPLICIT NONE
        REAL, INTENT(IN) :: x

        MySIN = SIN(DegreeToRadian(x))
     END FUNCTION MySIN

!   --------------------------------------------------------------------
!   FUNCTION MySIN():
!      This function takes a REAL argument in degree and computes its
!   cosine value. It does the computation by converting its argument to
!   radian and uses Fortran's cos().
!   --------------------------------------------------------------------

     REAL FUNCTION MyCOS(x)
        IMPLICIT NONE
        REAL, INTENT(IN) :: x

        MyCOS = COS(DegreeToRadian(x))
     END FUNCTION MyCOS

END MODULE MyTrigonometricFunctions
Click here to download this module. You also need a main program to test it. This mean program is identical to the
one used in a previous example. If you need it, click here to download a copy.

In this module, there are four PARAMETERs. Of these four, only PI is not listed as PRIVATE
and hence can be accessed from outside of this module. There are four internal functions,
MySIN(), MyCOS(), RadianToDegree() and DegreeToRadian(). The former two are listed as
PUBLIC and can be accessed from outside of this module. The latter two are listed as
PRIVATE and therefore cannot be accessed from outside of this module.

Note that if PI is also made PRIVATE, then the main program will have a mistake since it
displays the value of PI:

PROGRAM TrigonFunctTest
   USE MyTrigonometricFunctions

     IMPLICIT NONE
        ..........
     WRITE(*,*) 'Value of PI = ', PI                   ! PI cannot be used here
      ..........
END PROGRAM TrigonFunctTest

				
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
views:19
posted:11/20/2011
language:English
pages:53