Project 6 —Fortran 90

Document Sample
Project 6 —Fortran 90 Powered By Docstoc
					                                  Project 6 — Fortran 90
                                           David Arnold

                                          March 12, 2002


1    Introduction
It’s time to leave the old behind and discover the new. Fortran has undergone a number of revisions
since 1977 (hence the name, Fortran 77), and it’s time that we take advantage of the improvements
provided by Fortran 90.
    Unfortunately, the g77 compiler, provided by the Cygwin tools, will not compile a standard Fortran
90 program, so we will have to move our work to a new compiler. This is difficult, because there aren’t
many free Fortran 90 compilers available on the internet. Most of the compilers that handle modern
Fortran are prohibitively expensive. So, in order to save you some money, we are going to use the free
F compiler from Imagine.
    There a number of different ways that you can obtain the F compiler

    • You can download a free version of the F compiler at:

                                    ftp://ftp.swcp.com/pub/walt/F/

      You want to download the version appropriate for your operating system. For example, if you use
      the Windows operating system, then you want to download the file f win 020224.zip. Installation
      directions are availble in the file:

                               ftp://ftp.swcp.com/pub/walt/F/INSTALL

      You must have the GNU gcc compiler installed, but if you have already installed the Cygwin
      tools, the gcc compiler is already present on your system. When you’ve been using g77, the gcc
      compiler is called to do much of the work in compiling your Fortran 77 programs. If you do not
      have the Cygwin tools installed, then you need to download and install f win mingw 020224.zip.
      This is a DOS/Windows version of the GNU compiler tools.

    • You can purchase the F compiler on CD, along with the book, Programmer’s Guide to F , for
      $30. Information for this purchase can be found at:

                            http://www.fortran.com/imagine1/books.html

    • Finally, the software for the F compiler is also available on the CD we burned for you at the
      beginning of the semester. Be aware that several updates to the software have been made since
      we burned CD’s at the beginning of the semester. So, if you want the latest version of the F
      compiler, see me to burn a new CD, or download it from the site above.



                                                  1
1.1    Why F?
F is a programming language that serves as an introduction to programming for students in mathematics,
science, and engineering who have not programmed before. It is meant to replace such languages as
BASIC and Pascal, that in previous years, served to introduce students to the art of programming.
Typically, once a student finished an introductory programming course using either BASIC or Pascal,
then the student was prepared for the rigors of the far more difficult languages such as Fortran, C, and
C++.
    You might worry that you are learning a whole new language that is not related to Fortran, but let
me set your mind at ease, because this is not the case. Although the F programming language does
not contain all of the syntax or constructs available in modern Fortran 90 and 95 languages, F is a
proper subset of Fortran 90. This means that any program written in the F programming language will
compile on a standard Fortran 90 compiler.
    Therefore, as you study the F programming language, you are learning to program in Fortran
90. Think of the F programming language as an introduction to the Fortran 90 language. Thus, an
investment in the F programming language is an investment in your future programming needs in
mathematics, science, and engineering.

1.2    Installing The F Compiler
Instructions for installing the F compiler are available on the CD burned for you at the beginning of
class. However, we will take some time to explain the installation process in this article, and we will
compile a sample F program to test your installation.
    First, note that there is a new version of the F compiler available at the address:

                                  ftp://ftp.swcp.com/pub/walt/F

If you are a windows user, download the file f win 020224.zip.
    I encourage you to install the F compiler on your c: partition. If needed, move some things elsewhere
as F is a bit quirky when placed on other partititons. To start the installation, copy the program
unzip.exe from the CD to your c: drive, then copy the zip file f win 020224.zip to your c: drive as
well. Now, open a DOS box and change the directory to the c: drive with:

c:

Make sure that you are in the root directory of the c: drive with this command.

cd \

Now, type the following command to unzip the F distribution.

unzip f_win_020224.zip

This will create an F folder on your c: drive that contains the files for the F compiler to work. To
check this, open Windows Explorer, then browse to the c: drive and check out the contents of the F
folder. You should see subdirectories bin,doc, examples, lib, and tmp. The bin folder contains the
F compiler, which is what we will use to compile Fortran 90 programs. The examples folder contains
some excellent examples of Fortran 90 programs.
    You need to add the F\bin folder to your path. If you are using Windows 95 or 98, add the following
line to your autoexec.bat file, then reboot your system. See the Miktex Installation Page for a review
of this technique.


                                                   2
set PATH = %PATH%;C:\F\bin;
set TMPDIR=C:\F\tmp
In Windows 2000, add the above to the beginning of your path variable and add the TMPDIR environment
variable as shown in the Miktex Installation instructions on the CD burned for you at the beginning of
class.
    We now need to test the installation. Open a DOS box, then change directories with this command:
cd c:\F\examples
Next, compile the sample program seven 11.f95 with:
F -o seven_11 seven_11.f95
If all goes well, this will create a file called seven 11.exe in the examples directory. You can execute
this new file with the command:
seven_11
If all is well, the program should produce the output: "The percentage of rolls that are 7 or 11
is 23.30." Of course, since we are rolling dice, it is probable that you will get different numbers.

1.2.1      Getting F Documentation In Place
Open the Windows explorer (File Manager) and browse to c:/F/doc. There you will find documentation
files with extensions such as *.1 and *.3.
      • All files ending in *.1 need to be moved to the directory c:/cygwin/usr/man/man1.
      • All files ending in *.3 need to be moved to the directory c:/cygwin/usr/man/man3
Once you’ve moved all of these documentation files, try typing
$ man F
at the Cygwin prompt 1 . Note the uppercase F. That is, man f will not work.
    This command will open the manual pages for the F compiler, where you can read all about the
various options, bells, and whistles of the F compiler.

1.2.2      Adjusting .Emacs
There are a number of additions to make to your .emacs file to simplify the writing of Fortran 90 source
files in Emacs. Open c:/.emacs in Emacs and append the following lines.
;; Auto-fill

(add-hook ’fortran-mode-hook ’turn-on-auto-fill)
(add-hook ’f90-mode-hook ’turn-on-auto-fill)

;; Abbrev-mode

(add-hook ’fortran-mode-hook ’(lambda () (abbrev-mode 1)))
(add-hook ’f90-mode-hook ’(lambda () (abbrev-mode 1)))
These lines in your .emacs file will turn on auto-fill and abbreviation mode for both Fortran and Fortran
90 source files.
  1
      Note: In this command, the $ symbol is the Cygwin prompt.


                                                          3
1.2.3      Checking Errors with Emacs
By this time, you are probably used to typing C-x ‘ to go to the next error in your file. Unfortunately,
Emacs has not be set up to interact with the F-compiler, so we will have to make an adjustment if we
want this Emacs error checking behavior to continue.

   1. First, download the file:2

           http://online.redwoods.cc.ca.us/instruct/darnold/fortran/resources/compile.el

   2. Place the file compile.el in the directory c:/emacs/lisp/progmodes on your system.

   3. Now, open the file c:/emacs/lisp/progmodes/compile.el in Emacs. A new menu item will ap-
      pear named Emacs-Lisp. On this menu you will find a selection Byte compile this file. Click
      this selection with your mouse. This will compile compile.el, creating a file called compile.elc,
      which is stored as c:/emacs/lisp/progmodes/compile.elc.

   4. Close and restart Emacs.

From this point on, you can scroll through errors with C-x ‘. We will test this feature in the next
section.

1.2.4      Hello World Again
We now return to the first program we wrote in this class, but this time, we write the code in Fortran
90. Open a new file in Emacs as hello.f90. Note the extension is .f90 instead of .f. Emacs recognizes
this extension and automatically places itself in Fortran 90 mode. Enter these lines and save the file.

program hello
print *, ‘‘Hello, world!
end program hello

There are a couple of important notes to be made about the formatting of a Fortran 90 program.

       • First, note that Fortran 90 programs, such as hello.f90, are set in so-called “free format.” That
         is, no more reserving of the first six columns for comments, statement numbers, and line continu-
         ation characters. Of course, this is a relief, even though Emacs did that for us automatically.

       • Both upper and lower case characters are allowed, and variable names can be longer than six
         characters. Indeed, variable names can be up to 31 characters in length.

       • Although the official standard for Fortran 90 allows for the use of write statements, the F pro-
         gramming language does not allow the use of either write or format statements. However, the
         F programming language provides an easy-to-use print command that allows the output of the
         statement print *, ‘‘Hello, world!’’ to appear on your computer monitor. Soon, we will
         see how to replace the *, which is “free format,” with some format statements similar to those we
         wrote in Fortran 77.
   2
    As we go along, I learn more and more about setting up Emacs on computers. Thus, the download files that we burned
into your CD at the beginning of the semester are always in a state of flux. Files are added, updates are downloaded, etc.
You can always, if you wish, burn an additional CD to get the latest updates.




                                                           4
   • Note that the end statement of a Fortran 90 program is a little more complex than the end
     statement of a Fortran 77 program. Indeed, the keyword end must be followed by the keyword
     program and the same program name declared in the first line of the program. This construct is
     also required for subroutines and functions.
We’ve intentionally made an error in the source code of hello.f90, leaving off the double quotes at the
end of the string “Hello, world!” We’ve made this error on purpose, so that we can test the configuration
made in the section 1.2.3. We compile the file hello.f90 in the usual manner, only this time we call
the F compiler instead of g77.
F -o hello hello.f90
The compiler reports the following errors, the first of which is reached with C-x ‘.
F -o hello hello.f90
Error: hello.f90, line 1: Illegal character - octal value 15
***Malformed PROGRAM statement
Pressing C-x ‘ takes us to the first error, which unfortunately, is on line 1 of our program. Careful
scrutiny of the first line of our program reveals no visible problems. Once again, we suspect the CRLF
issue. This is indicated by the error message:
Error: hello.f90, line 1: Illegal character - octal value 15
The character with octal code (base eight) 15 is the DOS end-of-line character that must be changed
to a UNIX end-of-line character.

1.2.5   CRLF Again
We must return to our program source hello.f90 and save the file in UNIX format. If you recall, the
following sequence converts all end-of-line characters into UNIX format (C-x <ret> f unix <ret>).
That is, hit Control-x, followed by the Enter key, then press the letter f on the keyboard and type in
unix. Hit the Enter key to complete the command sequence.
    If all goes well, you will see the word (Unix) in the lower left corner of the status bar in Emacs. Be
sure to save the file. Now, try compiling again.
F -o hello hello.f90
We still get error messages and C-x ‘ takes us to the first error message.
Error: hello.f90, line 2: Unterminated character literal
       detected at ,@’Hello, world!
This is easily interpreted. We’ve forgotten the closing ” at the end of the string “Hello, world!” Fix
this error.
program hello print*, "Hello, world"
end program hello
Compile again. There should be no more errors.
cd c:/temp/ F -o hello hello.f90

Compilation finished at Wed Mar 13 17:04:03
Either in Cygwin, or in a DOS box, run the program as usual.
$ ./hello
 Hello, world

                                                    5
2     Introduction to Fortran 90
It’s now time to learn how to program in Fortran 90. Fortunately, much of Fortran 77 still works, but
there’s lots of new stuff to learn. Let’s begin.

2.1   Declaring Variables and Parameters
As in Fortran 77, Fortran 90 has integer, real, character, logical, and complex data types. There are
others, but we postpone talking about them for the moment. What’s different is the syntax used to
declare variables, but it’s not that intimidating and the syntax is easy to learn. For example, to declare
an integer variable ncount, we do the following.
integer :: ncount
Simple, the keyword integer, followed by two colons, then the variable name. The spaces before and
after each colon are not required, but they do help to make your source code easier to read. In a similar
manner, here’s how to declare a real variable.
real :: pi
Logical variables are also easily declared.
logical :: flag
    Fortran 90 allows you to declare a variable as a parameter constant and initialize the constant, all
in a single program statement. For example, suppose that we want to declare the real variable pi to be
a constant, with value 3.14159.
 real, parameter :: pi = 3.14159
Again, the spaces are not required, but they do help to make the source easier to read. In similar
fashion, you can declare other types of constants.
integer, parameter :: nCount = 100
logical, parameter :: flag = .true.
    Finally, here’s how you declare a character variable.
character (len = 80) :: text
Again, the spaces are not required, but they do help to make the source easier to read. Here is an
example of a character constant.
character (len = *), parameter :: &
    firstName = ‘‘David’’
Unlike Fortran 77, the length may be declared with an asterisk (even in the main program) indicating
that the length is to be determined from the value of the string. Furthermore, note that strings are
delimited with quotation marks instead of apostrophes.
    Finally, note the use of the line continuation character & in this character constant declaration.
This is the recommended use in most texts that I read, but Emacs has an annoying habit of using the
line continuation character twice when you break a line with C-c <ret>.
program hello
character (len = *), parameter :: &
     & msg = "Hello, world!"
print *, msg
end program hello

                                                    6
2.2     Control Structures
Some of the best improvements offered by Fortran 90 are its greatly simlified control structures.

2.2.1    If...Then...Else If...Else...End if
First, if...then...else if...else...end if works in precisely the same way as it does in Fortran
77.

program ifthenelse
  character (len = 20) :: firstName
  read *, firstName
  print *, "Your first name is: ", firstName
  if (firstName == "David") then
     print *, "Hello, ", firstName
  else if (firstName == "Adam") then
     print *, "Good bye, ", firstName
  else
     print *, "Hang up the phone, ", firstName
  end if
end program ifthenelse

Nothing new here, but now let’s take a look at do loops.

2.2.2    Do Loops
One important difference is the fact that a loop counter is not needed. This allows you to loop forever
or until some condition is met. For example, suppose that you want to find the largest power of 2 that
is smaller than 2000, but you are too lazy to do the math.3 A do loop with an exit statement will solve
this puzzle.

program pow2
  integer, parameter :: max = 2000
  integer :: i = 1, pow = 1
  do
     if (pow < max) then
         pow = pow*2
         i = i + 1
         print "(i10,i10)", i, pow
     else
         exit
     end if
  end do
end program pow2

There are several important points to be made with regard to this program.
   3
    Set 2n = 2000, then ln 2n = ln 2000 and n ln 2 = ln 2000. Thus, n = (ln 2000)/(ln 2) ≈ 10.96584. Thus, 210 = 1024 is
the highest power of 2 less than 2000.




                                                           7
   • First, note the do...end do structure. There is no continue statement in the F programming
     language. This also eliminates the need for a statement number at the end do position, which
     now replaces the former use of the continue statement.

   • Note the new exit command. When encountered in a do loop, this command exits the do loop,
     passing control to the program statement immediately following the closing end do statement.
     The exit command takes on added importance when one learns that the F programming language
     does not allows goto statements.

   • Finally, I could not resist putting in a little formatting of the output. Even though the F pro-
     gramming language has eliminated the write and format statements, as you can see, formatting
     the output is still a possibility in F, and it’s done in a natural way.
Here’s the output of the program.
          2           2
          3           4
          4           8
          5          16
          6          32
          7          64
          8         128
          9         256
         10         512
         11        1024
         12        2048
Subtle repositioning of the statements produces a “better” result.
program pow2
  integer, parameter :: max = 2000
  integer :: i = 1, pow = 1
  do
     if (pow<max) then
         print "(i10,i10)", i, pow
     else
         exit
     end if
     pow=pow*2
     i=i+1
  end do
end program pow2
And the output.
          1           1
          2           2
          3           4
          4           8
          5          16
          6          32

                                                  8
           7         64
           8        128
           9        256
          10        512
          11       1024
    Fortran 90’s cycle command allows you to return to the beginning of a do loop, forgoing execution
of the remaining statements in the loop.
program threemult
  integer :: i
  integer, parameter :: N = 20
  do i = 1, N
     if (modulo (i, 3) == 0) then
         cycle
     else
         print "(i5)", i
     end if
  end do
end program threemult
This program prints the numbers from 1 to 20, inclusive, skipping any multiples of three.
     1
     2
     4
     5
     7
     8
    10
    11
    13
    14
    16
    17
    19
    20
This result warrants some explanation.
    • The command modulo(i, 3) divides i by 3 and returns the remainder. Thus, if modulo(i, 3)
      equals zero, then the remainder when i is divided by 3 is zero, making i a multiple of three. When
      this happens, the cycle command returns control to the top of the loop, forgoing execution of
      remaining statements in the loop.


3    Program 6 — Integration
Often times, in calculus, you are presented with an integeral that has no closed form solution; that is,
a solution that can be expressed in terms of elementary functions such as polynomials, exponentials,
logarithms, trigonometric functions, and the like. This is a more common occurrence than one would

                                                   9
expect. Indeed, most of the time, you better have a numerical routine handy for doing your integration
because you are not going to be able to do the integral by hand.
    However, in this activity, we will choose an easily found definite integral for purposes of checking
the results of our program. Once you have confidence that your program is producing accurate results,
then you can use your program on integrals that have no closed form solution.
    So, with these thoughts in mind,
                                  2                   2
                                                1 4           1 4           15
                                      x3 dx =     x       =     (2 − 14 ) =    = 3.75.                (1)
                              1                 4     1       4             4
With this done, let’s discuss a method for finding a definite integral of an arbitrary function on an
interval [a, b].

3.1   The Left Endpoint Approximation
We first partition the interval [a, b] into n equal subintervals so that a = x0 < x1 < x2 < · · · < xn−1 <
xn = b. Next, use the left-endpoint of each resulting subinterval to calculate the height of a rectangle,
as shown in Figure 1. We sum the areas of the rectangles as follows,
                       y




                                                                                               f




                                                                                               x
                             x0        x1       x2                               xn   1   xn
                             a                                                             b

              Figure 1: Using the left endpoint to calculate the height of each rectangle.

                             S = f (x0 )∆x + f (x1 )∆x + · · · + f (xn−1 )∆x,                         (2)
where ∆x is the width of each rectangle and is calculated by dividing the length of the interval [a, b] by
n.
                                                    b−a
                                             ∆x =                                                      (3)
                                                      n
The width of each rectangle, ∆x, is a common factor, so we can write
                                  S = [f (x0 ) + f (x1 ) + · · · + f (xn−1 )]∆x.                      (4)
This will simplify the computation somewhat and ease the construction of the programming task. It is
S that you must compute to find the integral using the left endpoint method.

                                                              10
3.2   The Right Endpoint Approximation
We first partition the interval [a, b] into n equal subintervals so that a = x0 < x1 < x2 < · · · < xn−1 <
xn = b. Next use the right-endpoint of each resulting subinterval to calculate the height of a rectangle,
as shown in Figure 2. We sum the areas of the rectangles as follows,
                       y




                                                                                           f




                                                                                           x
                             x0    x1      x2                                xn   1   xn
                             a                                                         b

             Figure 2: Using the right endpoint to calculate the height of each rectangle.

                               S = f (x1 )∆x + f (x2 )∆x + · · · + f (xn )∆x,                         (5)
where ∆x is the width of each rectangle and is calculated by dividing the length of the interval [a, b] by
n.
                                                    b−a
                                             ∆x =                                                      (6)
                                                      n
The width of each rectangle, ∆x, is a common factor, so we can write
                                  S = [f (x1 ) + f (x2 ) + · · · + f (xn )]∆x.                        (7)
This will simplify the computation somewhat and ease the construction of the programming task. It is
S that we must compute to find the integral using the right endpoint approximation method.

3.3   The Midpoint Approximation
We first partition the interval [a, b] into n equal subintervals so that a = x0 < x1 < x2 < · · · < xn−1 <
xn = b. Next, use the midpoint of each resulting subinterval to calculate the height of a rectangle, as
shown in Figure 3. We sum the areas of the rectangles as follows,
                           x0 + x1        x1 + x2                xn−1 + xn
                  S=f              ∆x + f         ∆x + · · · + f           ∆x,                        (8)
                              2              2                       2
where ∆x is the width of each rectangle and is calculated by dividing the length of the interval [a, b] by
n.
                                                    b−a
                                             ∆x =                                                      (9)
                                                      n

                                                      11
                       y




                                                                                           f




                                                                                           x
                               x0   x1    x2                            xn   1   xn
                               a                                                  b

               Figure 3: Using the midpoint to calculate the height of each rectangle.

The width of each rectangle, ∆x, is a common factor, so we can write
                                x0 + x1        x1 + x2               xn−1 + xn
                    S= f                  +f             + ··· + f                    ∆x            (10)
                                   2              2                      2
This will simplify the computation somewhat and ease the construction of the programming task. It is
S that we must compute when evaluating the integral using the midpoint approximation.

3.4   Your Task
Write a program that will approximate the integral in equation (1). As this is your first Fortran 90
program, try not to use subroutines and functions, as I still have much to say about these Fortran 90
structures before you are ready to use them.
    Your program should do the following:
   • Query the user, asking which method the user wants to use:

       1. Left endpoint,
       2. Right endpoint,
       3. Midpoint, or
       4. all three results.

   • Once the user selects a routine, query the user as to the number of rectangles they wish to use.

   • Provide a nicely formatted output to the screen of the results of the user’s selection.

   • Finally, if you want to stretch the extra mile, gaining a little extra credit, add two more routines
     to your list.

       1. The Trapezoid Method.

                                                    12
                             http://online.redwoods.cc.ca.us/instruct/darnold/
                                    Math50B/CalculatorPrograms/trap.pdf
         2. Simpson’s Method.
                             http://online.redwoods.cc.ca.us/instruct/darnold/
                                    Math50B/CalculatorPrograms/simp.pdf


4     The Grading Rubric
The following rules will apply for this program, after which we will discuss and adjust the rubric during
class.

    1. (30 points) Will be awarded for adequate comments. Comments should include:

       (a) A description of the program’s purpose.
       (b) Name, date, version or revision number.
        (c) A complete dictionary of all variables and parameters used in the program.
       (d) Interprogram comments should proceed any code snippets explained by the comments. These
           should be adequately sprinkled throughout your code.

    2. (50 points) Will be awarded if the program works and does what it was asked to do.

    3. (10 points) Will be awarded for good program style. This includes good indentation practices,
       etc.

    4. (10 points) Will be awarded for creativity and extra effort. Did you just do the bare minimum?
       Or did you stretch and reach a little higher? Did you put something cute or clever into your
       program that nobody else seemed to think of?


5     Penalties
Each program that is assigned during the term will have a due date. On that date, the program must
be on the instructor’s desk before the start of class. Penalties will be assessed as follows.

    1. (10 points) There will be a 10 point deduction for any program that is handed in after the class
       has begun.

    2. (20 points) There will be a 20 point deduction per class period. That is, if you hand the program
       in one class period late, there is an automatic 20 point deduction. Two class periods warrants a 40
       point deduction, etc. To be clear, if the program is in the instructor’s hands before the beginning
       of the next class, that is a 20 point deduction. If the program is in the instructor’s hands before
       the start of the second class period past the due date, that is a 40 point deduction, etc.


6     Managing Files and Folders
Each of you has been given personal space on the sci-math server to store your work. Typically, this
space is mapped to the drive letter H. If you open the Windows Explorer (the file manager, not the
internet browser), you can see that the drive letter has been mapped to your login name.

                                                    13
   In this folder, create a new folder call FortranPrograms. Note that you must never use spaces in
filenames. In the Windows operating system, filenames are not case-sensitive, which is exactly opposite
what happens in Unix and Linux, where filenames are case-sensitive.
   In your FortranPrograms folder, create another folder called Program6. It is in this folder that
you are to place the source code and executables for this current project. When you receive your next
project, create a new folder called Program6 to hold that project, etc.
   If you work at home, I still want you to place copies of your work in the space reserved for you on
our system. Simply copy your home files onto a floppy disk and bring them with you to school. Use the
Windows Explorer to copy the files on your disk into the proper folder (H:/FortranPrograms/Program6).
   If everyone follows these simple rules, I can easily access your work from my office machine for
purposes of assigning a grade.


7    Caveat
On this project, if you stop by my office with hardcopy of your program before the due date of this
assignment, I will give a quick glance and critique of your source code. Somewhat like receiving a grade
on a draft before submitting your final draft for assessment.




                                                  14