During your 3 or 4 years as an undergraduate you will use computers for a number of tasks; experiment control, data collection, data analysis, image processing, computational physics, word processing. All of these require use of simple keyboard skills, text editing, simple operating system commands and software that other people have written. However as a physicist you will also want to programme computers yourself to solve problems and perform complicated tasks which you devise. This workshop is your first step towards becoming a competent computer programmer as well as a skilled computer user.
There are a bewildering number of types of computer, operating systems and programming languages and you will encounter several different systems as you progress through your degree course. However we have tried to design the first, second and third/fourth year computing workshops in a coherent way so that you gain a depth of experience in a system environment which is close to that used by professional physicists. So in this introductory workshop you will be using the same hardware and software as research workers in the University.
The computer service you will use is called Irix which is based on hardware produced by Silicon Graphics Inc. or more simply SGI. Their machines run a version of the UNIX operating system called IRIX. UNIX is a very widespread operating system which has been adopted by many manufacturers so when you are familiar with IRIX you will be able to use a large number of different machines. In the first and second year computing workshops you will learn Fortran. This is a high level computer language which is very suited to scientific computation and is in widespread use within the academic community around the world.
Access to the Irix system is via PC terminals that are situated in the Physics Department between the first and second year laboratories. These PCs run Windows 2000 operating from the CFS (Computer File Store).
In order to log on to Irix you will need a user name and password for both the CFS system and the Irix service.
As you progress though the workshop please keep notes about what you do. At the very least print out listings of your programs and stick them in your laboratory note book. You will find these useful later on in the course.
There are many good books on UNIX and Fortran but you don't need to purchase a book to do this workshop. In fact we recommend you do the workshop before considering buying a book on computing because after a little experience you will be in a better position to choose the right book for you. It is most likely that you want to buy a book on Fortran for use as a reference in the second and third/fourth year. Whichever book you choose make sure it is up to date and covers the 1990 standard. Recommended books are:
The scripts for this and other computing workshops and a guide to using Unix (Irix) can be found on the World Wide Web at URL:
http://www.star.le.ac.uk/compshop/
You must first log on to Windows 2000. This will bring up a typical PC window. A program called Exceed is used to connect to Irix. You get to Exceed from ``Start'' at the bottom left hand corner of the screen. Goto ``Program'' and then ``Exceed''. When Exceed starts it offers a choice of two machines harrier or owl. Clicking with the mouse on either of these will bring up a login screen for a machine. This includes the prompt:
login:
Type it your Irix username followed by return. Note that UNIX is case sensitive so you must be careful to check that the caps lock key is not engaged (usernames are in lower- case). Then enter your password:
Password:
Your password will not be echoed at the screen but if all is well you will be logged onto the system. The first time you log on you will enter an automatic process that forces you to change your password. Your new password must be at least 6 characters, 2 of which must be alphabetic and one of which must be numeric. Dictionary words and pure numbers are not allowed. Use something which is not easy to guess and don't forget it or write it down. Finally you will get a system prompt
$
This is the operating system soliciting for commands. Anything you type now instructs the computer to do something. Each command must be terminated by return. This is the signal to the machine to read the line of text you have typed and act upon it. To log off the Irix machine type exit (return).
$ exit
This will return you back to Windows 2000. Don't forget to log off the CFS system using ``Logoff'' available in the Start menu at the bottom left hand corner of the screen. Please use Logoff rather than Shutdown as the latter will reboot the machine and this takes some time.
Never leave a terminal logged on to an account unattended!!
All the information stored in the computer (usually on magntetic disk) is held in files. Many of these files are unformatted or binary files and can only be read by the computer. However some of the files are text files which can be listed at your terminal or printed out (for you to read) as well as read by the computer. A text editor is a program provided by the operating system which allows you to create and modify text files. There are a large number of editors available for UNIX systems. The vi editor (short for visual) is available on all UNIX systems but is designed for experienced programmers and can be rather confusing to beginners so we shall begin by using the eve editor which is a simple screen editor that is easy to use and more than adequate for Fortran programming.
0.5 1.28 1.0 4.39 2.0 19.6 3.0 45.5 5.0 123. 7.0 247. 9.0 408. 11.0 607. 13.0 840. 16.0 1271.
The following command will tell the system to start editing a file called flight.dat. Don't forget to type return at the end of the command.
$ eve flight.dat
The first time you type this the file won't exist so the editor will create it for you. The editor will automatically enter a screen editing mode. From now on everything you type will go directly into the file and text will appear on the screen, return will start a new line and delete will remove the character immediately to the left of the cursor. Go ahead and type in the table of numbers above. The exact format you use doesn't matter providing there are 2 columns separated by spaces. If you make a mistake you can use the arrow keys to the right of the main keyboard to move the cursor around the file so that you can add or delete characters as required. As you proceed you may wonder what you must type to get out of this process. Everything just seems to be entered into the file. The editor (and many other interactive programs) use special keys to control what they do. Many of these are above or to the right of the main keyboard (like the arrow keys) however the simplest way to get out of eve is to press the ctrl key and while holding this key down press z. Using control z in this way gets you out of the editor back to the system prompt. The editor will save the new edited version of the file overwriting any previous version. If you prefer you can use the PF10 function key to exit from the editor. This has the same effect as control z. If you want to abandon an editing session you can use the DO key to request a command prompt and then type quit. If your terminal has no DO key then use control b. Be careful, once you have typed quit the editing session is lost and can't be recovered.
When you have returned to the system prompt you can check that a new file has been created using the ls command which lists the names of your files.
$ ls
Don't panic if you accidentally exit from the editor and overwrite a good version of a file with a badly edited one. The eve editor saves the old version as the file name.bak (flight.dat.bak in this case) so that you can recover from any finger trouble by copying this file to overwrite the corrupted version.
$ cp flight.dat.bak flight.dat
You can list the file you have created with the editor by using the more command
$ more flight.dat
This copies the file to your terminal a page at a time. You must press the space bar to pass from one page to the next but for the above table a single page (screen) will be sufficient. If you spot an error in the file when you list it with more get back into the editor and correct it.
The eve editor is capable of many more sophisticated operations. When you are in the editor the help key will provide a summary of what other special keys do. If there is no help key use control b to get a command prompt and then type help.
You can obtain a hardcopy listing of any text file by using the lp command.
$ lp -dA4 <flight.dat
I bet you wouldn't have guessed that. UNIX is notorious for being obscure, easy when you know the commands, impossible if you don't. The -dA4 is an option which sends the file to the A4 print queue. By default lp takes its input from the standard input which is your terminal, so <filename changes the standard input to the file you want to print (< is a redirection operator in UNIX speak). In order to get your hardcopy you must obtain a Print Copycard from the demonstrator and use the Print Station in the terminal room between the 1st and 2nd year Laboratories. Don't print text files until you are certain they are correct.
Now you can use the editor you are ready to create a Fortran source file. This is a text file which contains instructions for the Fortran 90 compiler and by convention has a name which ends in .f90 on UNIX systems.
$ eve first.f90
Enter the editor using the above command and type in the following simple Fortran program.
PROGRAM FIRST
! This is a simple Fortran program
IMPLICIT NONE
REAL :: R1,R2,SUM
PRINT *, "Type 2 numbers"
READ *, R1,R2
SUM=R1+R2
PRINT *, "The sum is ",SUM
END PROGRAM FIRST
Fortran 90 is written in free form and so the details of the layout are upto you but be careful to type the correct characters. The text enclosed in double quotes, "Type 2 numbers" and "The sum is " are literal strings which are read just as a string of characters by the compiler without any alteration or interpretation. When you have completed typing the program in, exit from the editor and check that you have a file called first.f 90 using the ls command and list it on your terminal using more.
$ ls $ more first.f90
If you're happy that the source file looks correct then you can try using the compiler.
$ f90 first.f90 -o first
The compiler reads the source text file and translates it into a binary form called the object code. Because you have written a complete program the f90 command then takes this object code and links it together with all the system object code that is needed to make the program run. This produces an executable module which is written to the file first that is specified by the -o option on the command line. If your source file is incorrect then the compilation will fail and error messages will tell you what the problem is. In this case you must re-edit the source file to correct the mistake and try the f90 command again. When things are alright a file first will be produced. You can then run the program by just typing the file name first.
$ first
The program should ask you to type 2 numbers. If you type them in (separated by a space, comma or return) the program will list the sum and stop. If the result is wrong then the syntax of the program is correct (compiles ok) but the logic is not (gives the wrong answer!).
Programs are started by typing the name of the executable file. Once they are running they are in control but here are various control keys you can use to control what goes on.
ctrl C to interrupt the process and return to the $ prompt ctrl Z push the process into the background (stopped) ctrl D indicate EOF (end of file) at the terminal
In each case you must hold down the ctrl key and then press the letter key. If you accidentally use control Z you can bring the process back to the foreground using the command fg. Processes that are in the background can be listed using the command jobs.
$ jobs $ fg
You should use control C if your program is stuck in an infinite loop or refuses to stop for some other reason.
Control keys can also be used to recall and edit previous command lines:
ctrl P recall previous command lines ctrl N move to next command line (opposite of ctrl P) ctrl B move cursor back along the current command line ctrl F move cursor forward along the current command line
In the introduction you typed in a very short Fortran 90 program consisting of a number of lines or statements. Apart from the line starting with ! each statement is an instruction to the compiler and it is usual to start each statement on a new line. The program declared and used variables, R1, R2 and SUM. These variables are the way Fortran (and in fact all high level languages) refer to locations in the computer's memory. The program FIRST also used READ and PRINT statements. These are the way Fortran performs input and output (usually refered to as I/O) with the outside world. The sections below introduce you to the common forms of Fortran statements. As you read through the sections you will be requested to type in and test the example programs given to try out the statements as they are introduced. There is a lot of information in these sections. Take your time and read it carefully.
Fortran 90 compilers recognize the following characters :
26 letters of the alphabet A-Z uppercase only (Fortran compilers ignore case) 10 decimal digits 0-9 4 arithmetic symbols + - * / (note ** for raise to the power) 18 other symbols ( ) . = , ' $ : space ! " % & ; < > ?
All compilers will accept other characters in comment fields but these are ignored in any compilation. It is good style to use uppercase for Fortran statements and the usual mixture of upper and lower case for the comment field.
Fortran 90 accepts statements in free form in which the detailed layout is unimportant. More than one statement can occur on a single line using a semi-colon to separate them. However the syntax (spelling, puncuation and order) of the statements is most important.
If ! (exclamation mark) is found then the rest of the line is taken to be a comment and is ignored by the compiler (look at the first program above). Sometimes a statement is too long to fit on one line. When this occurs an & (ampersand) can be put at the end of the line to indicate that the next line should be joined on as a continuation of the current line.
Many types of numbers and other entities are manipulated by computers, integers (whole numbers), floating point (real) numbers consisting of a mantissa and exponent, (the mantissa holds the significant digits and the exponent holds the scaling factor or power of 10) encoded character strings, unsigned bytes etc.. There are three forms of quantity manipulated by Fortran, literals, variables and parameters. Literals are text strings which directly represent the quantities concerned. For example typical literals are:-
Integer constants 1 411 -34 +18 Real constants 1.35 44.33 -0.0123 1.3E-3 Character strings "This is a character string"
In order to include the double quotes character in a literal text string you can use single quotes to delimit the string. Variables are quatities that can be changed during the execution of a program. Parameters are fixed quantites.
Variables and parameters have names. In Fortran we MUST make it clear to the compiler what type of quantity is associated with each name. This is done by declaration statements which are placed BEFORE any executable statements at the top of the program.
In Fortran the variable and parameter names must conform to the following rules:
up to 31 characters, alphabetic, underscore or numeric the first character must be alphabetic
So in the first program above the variables R1, R2 and SUM are real (floating point) numbers. When the computer writes them out it includes a decimal point. The declaration statements should start with an IMPLICIT NONE statement. If it is used ALL variables etc. must be explicitly declared at the start of your program. It is often important that you use the correct type of variable for a given purpose. Integers are usually used for counters, subscripts etc. whereas real numbers normally represent measured or continuous quantities.
It is often useful to refer to a constant (for example p ) using a name. Fortran uses parameter names to do this.
So typical declarations look like:
IMPLICIT NONE
REAL :: RR, ZZ
INTEGER :: I,IFLAG
CHARACTER(LEN=40) :: BLURB
INTEGER, PARAMETER :: NMAX=100
The specification of the declared type occurs to the left of the :: and the names occur to the right.
Let's look at the first program again.
PROGRAM FIRST
! This is a simple Fortran program
IMPLICIT NONE
REAL :: R1,R2,SUM
PRINT *, "Type 2 numbers"
READ *, R1,R2
SUM=R1+R2
PRINT *, "The sum is ",SUM
END PROGRAM FIRST
The PROGRAM statement declares the start of a new program section. Just 1 such section is required for each linked executable program. The name of the program is up to you. The next line is a comment line starting with ! which is ignored by the compiler. The next two lines declare the real variables.
The PRINT transfers data values from the program to the outside and the READ transfers data values from the outside into the program. The *, specifies that the input/output is list directed, that is the format of the text outside the program is determined by the list of quantities to the right of the comma. On READ the computer will try its best to read what is required whatever the format outside (from your terminal) and on PRINT the computer will choose the format it thinks best. The list of things to the right of the *, is the I/O list. It tells the computer which memory locations to READ and PRINT to/from. Note that elements within an I/O list are separated by commas. You can test out the READ and PRINT by changing the magnitude and type of the number you type in when first is run.
The END statement tells the compiler where the PROGRAM finishes. Every PROGRAM statement should be matched by an END statement.
Each READ or PRINT will normally deal with at least 1 complete line of text. If a READ reaches the end of a line before all the I/O list items have been satisfied then the READ will proceed to the next line. If there are more values on the line than required then the excess values will be ignored. By default every READ starts a new line. Similarly every PRINT will start a new line.
In the program first SUM is assigned the value of R1+R2. Such assigments take the form:
X=A+B ! set variable X equal to A+B
J=J-1 ! decrement variable J by 1
Y=X**2 ! set variable Y to X squared
Note that there is an implied flow from right to left and they are NOT the same as algebraic equations. Therefore incrementing or decrementing a variable is allowed. The operators for arithmetic etc. are, in order of precedence:
( ) parentheses to dictate the order of operations ** exponentiation ( raising to a power) * or / multiplication or division + or - addition or subtraction
You can mix numerical types within assignments and arithmetic statements but you MUST be careful. If you assign a REAL result to an INTEGER then truncation or rounding down will occur. Truncation will also occur if you perform INTEGER division.
I=R+35.3 ! The REAL result is truncated before assignment to INTEGER A=I/L ! The result is truncated to an INTEGER and assigned to REAL R2=R1+II ! The result of mixed arithmetic is REAL
Type in the following program which can read in the data file flight.dat you created when learning to use the text editor.
$ eve second.f90
PROGRAM SECOND
! Declare variables
IMPLICIT NONE
REAL :: T,H
INTEGER :: J
! Open a data file
OPEN(UNIT=1,FILE="flight.dat",STATUS="OLD")
! Read data from file and print it out
DO J=1,10
READ(UNIT=1,FMT=*) T,H
PRINT *, T,H
END DO
CLOSE(UNIT=1)
END PROGRAM SECOND
Compile and run the program. You should find that the values from the file flight.dat that you typed in earlier are listed at your terminal.
$ f90 second.f90 -o second $ second
The OPEN statement makes a connection between the file specified by FILE="flight.dat" and the unit number (sometimes called a channel number) specified by UNIT=1. The unit number is then used to refer to the file in the program, for example READ(UNIT=1,FMT=*). Now the READ will be from the file rather than your terminal. The STATUS="OLD" in the OPEN statement tells the computer to find a file that already exists. You must use STATUS="NEW" to create a new file. In that case you must, of course, WRITE to the appropriate channel rather than READ because the file will be empty! The CLOSE statement breaks the connection between the file and the program and releases the unit number for subsequent use. You should always CLOSE units when they are no longer required. You can choose any number for your channels (within reason!) but should avoid 5 and 6 since these are usually the default for the terminal input and output.
The program second also contains a DO loop. This is a repetition structure. The DO J=1,10 specifies a loop counter J and a range 1,10. Each time the DO is encountered J is incremented by 1 and when END DO is reached control is passed back up to the DO statement. So the loop is executed 10 times for values of J=1,2,3,4,5,6,7,8,9,10. If an increment (or decrement) other than 1 is required this can be specified as a third argument after the = sign.
DO K=1,13,2
other statements
END DO
If the loop is incrementing then control will pass to after the ENDDO when the loop variable (K above) is greater than the second argument. If the loop is decrementing then the loop will terminate if the variable is less than the second argument. If this condition is true at the start of the first pass, the loop will be skipped altogether.
Notice that the statements within the DO loop have been indented further along the line compared with the DO and END DO statements. This is not mandatory but does make the code much easier to read and understand. It is recommended that you indent all your structures in this way.
Subscripted or array variables can be declared by including a dimension specification in the declaration statement.
REAL, DIMENSION(100) :: AR
The name (AR in the above example) then refers to the collection of all the elements. Single elements are then referenced by using INTEGER subscripts enclosed in parentheses.
! Array elements are referenced using integer subscripts
A(1)=5.5
TABLE(1,2)=5
I=2
A(I)=13.0
J=4
R=A(J+I)
When declaring a character variable the length or number of characters must be specified.
CHARACTER (LEN=40) :: BLURB
We can improve our second program by declaring a character variable to hold the file name and declaring arrays to read the data into.
PROGRAM SECOND
! Declare variables (before executable code)
IMPLICIT NONE
CHARACTER (LEN=40) :: FNAME
REAL, DIMENSION(10) :: T,H
INTEGER J
! Get file name from terminal
PRINT *, "Type data file name"
READ *, FNAME
! Open file
OPEN(UNIT=1,FILE=FNAME,STATUS="OLD")
! Read data from file into arrays
DO J=1,10
READ(UNIT=1,FMT=*) T(J),H(J)
ENDDO
CLOSE(UNIT=1)
! List the data arrays
DO J=1,10
PRINT *, T(J),H(J)
END DO
END PROGRAM SECOND
Now the program will still work even if you change the data file name. After the first DO loop all the data from the file are held in the arrays T and H. They can be listed even after the data file has been CLOSEd.
Try modifying your version of program second to use arrays as indicated above.
We can enhance the above program still further so that it will work on data files with different numbers of lines (T,H data pairs). Copy your second.f90 source file to a new third.f90 file and try editing it to produce the following program.
$ cp second.f90 third.f90
$ eve third.f90
PROGRAM THIRD
! Declare variables (before executable code)
IMPLICIT NONE
CHARACTER (LEN=40) :: FNAME
INTEGER, PARAMETER :: NMAX=900
REAL , DIMENSION(NMAX) :: T,H
REAL :: TT,HH
INTEGER :: IERR,J,K
! Get file name from user
PRINT *, "Type data file name"
READ *, FNAME
! Open file
OPEN(UNIT=1,FILE=FNAME,STATUS="OLD")
! Read data from file into arrays
IERR=0; J=0; TT=0.0; HH=0.0
DO
READ(UNIT=1,FMT=*,IOSTAT=IERR) TT,HH
IF(IERR==0) THEN
J=J+1
IF(J>NMAX) THEN
PRINT *, "too many data points in file ",&
FNAME
J=J-1
EXIT
ELSE
T(J)=TT
H(J)=HH
ENDIF
ELSE
EXIT
ENDIF
END DO
CLOSE(UNIT=1)
PRINT *, J, " data points read from ",FNAME
! List the data arrays
DO K=1,J
PRINT *, T(K),H(K)
END DO
END PROGRAM THIRD
In order to cope with varying numbers of data points the arrays T() and H() must be declared large enough. We have used a PARAMETER declaration statement to declare the value NMAX. This is then used throughout the program. If we want to change the maximum size we just have to change this parameter definition. The definite DO loop has now been replaced by an indefinite DO loop controlled by the IF, ELSE, ENDIF structure. The (IERR==0) is called a conditional expression and the == is a relational operator. There are two forms for each of the relational operators for compatibility with older versions of Fortran, The complete set is :
== or .EQ. equals
/= or .NE. not equals
> or .GT. greater than
< or .LT. less than
>= or .GE. greater or equal
<= or .LE. less or equal
These operators compare 2 quantities and return a logical value, true or false, yes or no. The indefinite DO will execute another pass each time END DO is reached until IERR is not equal to zero. The integer IERR is set by the READ statement which specifies IOSTAT=IERR. This tells the computer to return the status of the READ to the INTEGER IERR. This status value is 0 if the READ is executed correctly. If some error occurs then the status has a positive value (related to the type of error). If the READ hits the end of the file (EOF) status is returned as -1. Within the loop we have used 2 decision structures. These take the form of IF, ELSE, ELSEIF, ENDIF statements. The IF and ELSEIF are always accompanied by conditional expressions. The lines which will be executed under each condition have been indented in so you can see the structure of the logic. Each IF must always be terminated by ENDIF and each DO must be terminated by ENDDO. When IERR is not equal to zero then the EXIT statement is reached. This forces a jump to the first statement after the END DO and therefore terminates the indefinite loop. Notice that an EXIT is also reached if there are too many data points to fit into the declared arrays.
It is possible to combine the logical results arising from the use of relational operators using logical operators.
IF(E1>E2.AND.E3<E4) THEN
statements
ENDIF
Note that the relational operators take precedence over logical operators. Other logical operators are:
.NOT.
.AND.
.OR.
.EQV.
.NEQV.
They are all delimited by dots just like the older forms of the relational operators.
Note that the order of the IF, ELSEIF , ELSE, ENIDF, DO and ENDDO statements is important. Such structures must not cross but they can be nested one inside the other as indicated by the layout of the indentation. This form of layout is not required by the compiler but makes the code much more readable for the programmer.
If repetition is required within a single statement, for example using a single READ to read many elements into an array, then it is often easiest to use a special syntax called an implied loop.
REAL , DIMENSION(10) :: A
READ(UNIT=1,FMT=*) (A(K),K=1,10)
Such a loop is always enclosed by parentheses and the commas must be included as shown.
You MUST NOT assume that variables will have a value when a program starts up although on many machines they may be set to zero by default. This is why in third above IERR and J are both assigned before the DO loop. If you require initial values to be set before the program starts this can be done as part of the declaration. Actually we have already seen this in the PARAMETER declaration of NMAX in the third program above.
REAL :: A=1.45
REAL, DIMENSION(2) :: B=(/ 1.23,3.45 /)
REAL, DIMENSION(50) :: C=(/ (0,I=1,50) /)
REAL, DIMENSION(40) :: D=0.0
Implied loops can be used to save space. In the last example above the array name D refers to all the elements of that array so the assignment is carried out on an element by element basis. Of course you can always use assignments to achieve the same thing although these are executed at run time rather than during compilation. Constants like p or the number of seconds in a day which may occur many times within a program should be assigned using a PARAMETER declaration. Never type such constants explicitly into the code as literal constants. Sooner or later you will make a mistake and the constant will change its value without you noticing. For example :
REAL, PARAMETER :: DAYSEC=86400.0
There are a number of mathematical functions which are built into Fortran 90. The arguments of these functions are enclosed in parentheses and multiple arguments must be separated by commas. They are used within expressions as follows:
Y=SQRT(X)
Z=MOD(Y,2.0)
Other common functions available are:
ABS absolute value
SQRT square root
SIN sine (argument radians)
COS cosine (argument radians)
TAN tangent (argument radians)
ASIN arcsine (returns radians)
ACOS arccosine (returns radians)
ATAN arctangent (returns radians)
ALOG10 logarithm base 10
EXP raise e to power
INT truncate REAL to INTEGER
NINT fix REAL to nearest INTEGER
REAL float INTEGER to a REAL
MOD remainder of A1/A2
Try modifying your second program to list the square root of the height of the rocket as a function of time from the file flight.dat What do you notice about the motion of the rocket? Can you predict the height after 20 seconds?
If an array name is referenced without indices indicating a single element then the whole array is implied. If this occurs on the left of an assignment then an implicit DO loop will be performed and all the elements will be assigned in turn. In order for this to happen the left and right hand side of the assignment must be conformable. This will be true if both are arrays with the same dimensions. Scalars are always conformable with any array. The following statements use this powerful convention which avoids unnecessary use of explicit DO loops.
REAL, DIMENSION(100) :: A,B,C
A=1.4 ! set all elements of array A to a scalar
B=A+4.0 ! add 4.0 to all elements of A
C=A*B ! multiply elements of A by elements of B
D=SIN(C) ! take the SIN of each element of C
As you can see the implicit loops also apply to the intrinsic functions. If the function requires more than 1 argument then the arguments must be conformable. In general if a statement seems sensible and there is no ambiguity then such whole array operations will work.
If you have worked through the text above you should have typed in 3 programs first, second and third which introduce elementary Fortran syntax and structures. You are now ready to write your own programs but please take note of the following advice about programming:-
An algorithm is a prescribed set of instructions for solving a problem. The following problems are designed to give you practice in writing simple algorithms. You should first try to work out how you would solve the problem without programming the computer using, for example, a calculator. Work out what the steps of the algorithm are and then decide how to translate these into Fortran 90. Endeavour to make the programs as neat, efficient and understandable as possible.
The first three problems involve just natural numbers, positive integers. The next three problems involve floating point (REAL) numbers. The last four problems involve analysis of the rocket flight data.
When writing algorithms there is usually a tradeoff between the time taken and the space (memory) required. This will be illustrated in your answers to problems 1), 2) and 3). The program required for problem 1) is efficient on space (no arrays are declared) but is slow. The programs 2) and 3) are faster but require a large array to store the prime numbers found.
The floating point numbers in any computer are only approximate. They are represented by two integers, the mantissa and the exponent. The mantissa holds the significant digits and the exponent holds the scale (power) usually in base 2. The mantissa is always normalised so that the most significant bit is set and the exponent is adjusted accordingly. This insures maximum precision within the space available.
Do this by the brute force method without saving the prime numbers as you find them. First get your program to request the upper limit (NMAX say). Then use a definite loop to work through each candidate number in turn. Inside this definite loop you will need an indefinite loop to look for factors. The function MOD(I1,I2) can be used to check if I2 is a factor of I1. If it is then this function will be zero.
The chosen number N should be declared as a PARAMETER in the program. You should declare an integer array to store the primes as you find them. For example to find the first 1000 primes your declarations could include:
INTEGER, PARAMETER :: N=1000 INTEGER, DIMENSION(N) :: IP
Since you only need to check for prime factors, the accumulating list of primes (stored in the array IP) can be used to test for factors of the next candidate. Don't list the primes until they have all been found. You will probably find it convenient to seed the list of primes with 2 to get things started (i.e. set IP(1)=2). Note that, apart from 2, all primes are odd so that every even candidate value can be skipped.
First get the program to request the target integer. Then search for prime numbers and when you find one check to see if it is a factor of the target integer. Note that the default size of integers stored by Fortran 90 is 4 bytes (32 bits) so the maximum positive integer which can be stored is 231-1 = 2,147,483,647.
I suggest that to start with you limit your algorithm to considering the first 3000 prime numbers. Use your program to find the prime factors of the integers:
630,210,767 1,646,913,833 2,108,801,987
You will now appreciate that a problem which is easy to state in words is not so easy to solve by writing an algorithm. The solution to problem 3 above is important. It is impossible to speed up the search if the target is large and if there are only a few prime factors. If the number becomes too large then it is impossible to solve in a sensible time even with the fastest machines available. Such large numbers can be used for making very secure codes. Problem 3 must be solved to break such codes.
(1+x+x2+x3+...xN)
N->infinity, x<1.0
Compare your answer with the analytical result
SUM=1/(1-x).
The number of terms required depends on the accuracy you want. Each successive term in the series is smaller than the previous term so check for convergence by calculating the fractional change afforded by each term. How many terms are required to yield a fractional error of order 10-5 if X=0.5?
(x-x3/3!+x5/5!-x7/7!...)
Compare your answer with the intrinsic Fortran 90 function SIN(X).
In this case you should check for an absolute accuracy of 10-5 in the result.
The problems above demonstate that the precision of floating point calculations is most important. Fortran 90 allows you to choose the precision of the floating point numbers.
INTEGER, PARAMETER :: PEND_PREC=SELECTED_REAL_KIND(P=10,R=30) REAL(KIND=PEND_PREC) :: SUM,TERM,DIFF
The function SELECTED_REAL_KIND(P=10,R=30) returns an integer which is used to specify the precision. P=10 requests at least 10 decimal digits in the mantissa and R=30 requests an exponent range of 1030. The REAL(KIND=PEND_PREC) then uses the integer returned to specify the precision of the declared variables. Try modifying your program to use this higher precision. Any REAL constants you use should be declared using D for the exponent.
1.0 a real floating point constant with single precision 1.0E0 a real floating point constant with single precision 1.0D0 a real floating point constant with double precision
! Loop to calculate approximate acceleration at successive times
! using finite differences. Also find mean acceleration.
AA=0.0
TT=0.0
DO J=2,ND-1
VL=(H(J)-H(J-1))/(T(J)-T(J-1))
TL=(T(J)+T(J-1))*0.5
VH=(H(J+1)-H(J))/(T(J+1)-T(J))
TH=(T(J+1)+T(J))*0.5
A=(VH-VL)/(TH-TL)
PRINT *,T(J),A
AA=AA+A*(TH-TL)
TT=TT+(TH-TL)
END DO
PRINT *,"mean accel. from finite differences ",AA/TT
y = mx + c
Given values of y corresponding to values of x we want to find the gradient (or slope) m and the intercept c. If the y values contain errors we want the best estimate of m and c. This can be done by least squares fitting (also called linear regression) a method fully explained in the laboratory script 'Error estimation in first year laboratory experiments'. You will need lines like the following.
! Perform linear regression on sqrt(h) vs. time
SUMX=0.0
SUMY=0.0
SUMXX=0.0
SUMXY=0.0
DO J=1,ND
SH=SQRT(H(J))
SUMY=SUMY+SH
SUMX=SUMX+T(J)
SUMXX=SUMXX+T(J)*T(J)
SUMXY=SUMXY+T(J)*SH
END DO
RN=REAL(ND)
DEL=RN*SUMXX-SUMX**2
GRAD=(RN*SUMXY-SUMX*SUMY)/DEL
C=(SUMXX*SUMY-SUMX*SUMXY)/DEL
What is the best estimate of the acceleration of the rocket? How does this compare with the results from problem 7? What is the best estimate of the launch time of the rocket (i.e. t for h=0)?
AC=2.0*GRAD**2
SDIFF=0.0
DO J=1,ND
PREDH=0.5*AC*T(J)**2
DIFF=H(J)-PREDH
PRINT *, T(J),DIFF
SDIFF=SDIFF+DIFF**2
ENDDO
RMS=SQRT(SDIFF/RN)
PRINT *,"rms height error from linear regression ",RMS
/disk/s/zrw/under/flight
is a larger version of the file you have typed in. Try out your software on this data file. If you have welded the file name into your code or only allowed for 10 data pairs you'll have to modify your code. The file contains 105 data points. This is a good test of how general and adaptable your programs are. Is the acceleration constant? Try copying the data file to your directory using the cp command and splitting the file into sections using the editor or the csplit command so that you can analyse each section separately. If you need help with this ask a demonstrator. What is happening to the rocket?
The Fortran described above is more than sufficient to tackle the problems 1 to 10 but it does not reveal the full power of the language. In particular we have not considered the modular structure possible in Fortran or the complications of formatted and unformatted I/O. These and other topics must wait until next year. However the remaining sections below give you a foretaste of what is to come and may be of interest if you want to write more complicated programs in your own time.
In Fortran the task the computer is to perform can be divided up into separate modules or program sections. Each such section is compiled completely independently although several sections can reside in a single source file. So far we have seen just one type of section, the PROGRAM. This is entered from the operating system when the image is started and returns control to the operating system when it stops. There are 2 other forms of program section, the SUBROUTINE and the FUNCTION. These can only be CALLed or entered from within a PROGRAM or another SUBROUTINE or FUNCTION. When they finish control is returned to the point at which they were CALLed. Both SUBROUTINEs and FUNCTIONS can have an argument list which allows values from one program section to be transfered from the calling routine to the called routine at the call time and returned back to the calling routine on the return. The declaration of a SUBROUTINE or FUNCTION to the compiler takes the form:
SUBROUTINE DOIT(A,B)
IMPLICIT NONE
REAL :: A,B ! declare dummy arguments
internal declarations
executable instructions
END SUBROUTINE DOIT
FUNCTION ZAP(C)
IMPLICIT NONE
REAL :: ZAP,C !declare function and dummy arguments
internal declarations
executable instructions
ZAP=expression !assign value to function before return
END FUNCTION ZAP
In such code the arguments in parentheses (A,B) or (C) are said to be dummy arguments. This is because the compiler knows that they actually refer to memory associated with the calling routine. Since this is unknown at the time of compilation they are not allocated memory within the current routine. However the type of such arguments must be made clear to the compiler. You do this using explicit declarations. Any variables which do not appear within the argument list are called local variables. They have no connection with variables elsewhere even if the names are the same. All program sections must finish with an END statement and in the case of SUBROUTINEs and FUNCTIONs there is an implicit RETURN just before the END.
Subroutines are invoked using a CALL. In the following example DOIT calculates a quadratic expression in X amd returns it as Y.
REAL :: X,Y
X=4.0
CALL DOIT(X,Y)
PRINT *, Y
.
.
SUBROUTINE DOIT(A,B)
IMPLICIT NONE
REAL :: A,B
B=4.0*A*A+3.0*A+10.0
END SUBROUTINE DOIT
The order and type of the argument list in a CALL must agree with the declared argument list for the routine. It is customary to declare input (or imported) arguments to the left of the list and output (or exported) arguments to the right with those (if any) arguments which are imported, modified and then exported in the middle.
Function subprograms declared in the above manner are used in the same way as the instrinsic functions. Note that there is a data type associated with the value returned by a FUNCTION and the value is returned by using an assignment to the FUNCTIONs name within the FUNCTION definition.
Care must be taken when passing arrays as arguments. The safest way is to include the dimensions as arguments.
SUBROUTINE THINGY(N,A,B)
IMPLICIT NONE
INTEGER :: N
REAL, DIMENSION(N) :: A
REAL :: B
Using this technique ensures that the routine knows the size of the array concerned and the dimensions are available as variables within the routine, for use in DO loops etc..
The exact format of a READ, PRINT or WRITE can be specified using a FORMAT statement. This can be referenced by a line number as the argument of the READ, PRINT or WRITE in place of the * or FMT=*.
PRINT 10, II,RR
10 FORMAT(5X,I3,1X,F10.3)
Alternatively the format can be embedded within the I/O statement as a character string:
PRINT "(5X,I3,1X,F10.3)", II,RR
The format is specified by a series of descriptors which are enclosed in parentheses after the word FORMAT. These descriptors are separated by commas. There are a large number of descriptors available but the 5 most common are given below. In each case the n specifies the total number of characters in the field of the descriptor.
nX n spaces
Fn.m REAL number, fixed point, m digits after point
En.m REAL number with exponent, m digits after the point
In INTEGER
An CHARACTER string
The F, E, I and A formats tell the computer how to READ,PRINT or WRITE the specified types. The complete FORMAT is scanned from left to right. Each descriptor requiring a type will expect a corresponding item in the I/O list of the READ, PRINT or WRITE statement with the correct type. The total width of the floating point number fields must include the decimal point, any negative sign etc.. All characters required to write out a number must be included. If the number is too small to fill the field of n then the result will be right justified padding out with spaces to the left.
The first character generated by the format on WRITE is used for carriage control and pagination unless you include extra keywords in the OPEN statement for the UNIT number.
space start a new line
0 advance 2 lines
1 start a new page
+ overwrite current line
The print statement automatically includes a space when the default output device is a printer but when using a READ statement to a printer (or a device configured like a printer) then the above control characters should be included at the beginning of the format:
98 FORMAT(1X,F10.5,1X,F10.5) 99 FORMAT(2(1X,F10.5))
26 letters of the alphabet A-Z uppercase only (Fortran compilers ignore case) 10 decimal digits 0-9 4 arithmetic symbols + - * / (note ** for raise to the power) 18 other symbols ( ) . = , ' $ : space ! " % & ; < > ?
; semi-colon is used to separate statements on a line ! exclamation indicates the start of comment (rest of line ignored) & ampersand at end of line indicates continuation on next line
Integer constants 1 411 -34 +18 Real constants 1.35 44.33 -0.0123 1.3E-3 Character strings "This is a character string"
variable names up to 31 characters alphabetic, underscore or numeric first character must be alphabetic without IMPLICIT NONE or some other IMPLICIT declaration then: if first character I,J,K,L,M,N then type INTEGER otherwise type REAL
All declarations within a code module must appear before any executable statements.
IMPLICIT NONE ! no name convension for type
INTEGER :: I,IFLAG ! declare integer variables
INTEGER, PARAMETER :: INUM=100 ! declare integer parameter (constant)
INTEGER, PARAMETER :: D_P=SELECTED_REAL_KIND(P=10,R=30)
! define precision parameter
REAL :: RR, ZZ ! declare real variables
REAL(D_P) :: DD,FF ! real variables with precision
REAL, DIMENSION(3) :: POS ! array dimension
REAL :: TIM(3) ! array dimension declaration
LOGICAL :: FLAG ! logical variable (true or false)
CHARACTER(LEN=40) :: BLURB ! character variable
CHARACTER(LEN=10), DIMENSION(3) :: TARR
! character array of 3 elements
! each element 10 characters
REAL, ALLOCATABLE, DIMENSION(:) :: X,Y
! dynamically allocatable arrays
REAL :: A=1.45 ! declare a real with initial value
REAL, DIMENSION(2) :: B=(/ 1.23,3.45 /)
! initial array element values
REAL, DIMENSION(50) :: C=(/ (0,I=1,50) /)
REAL, DIMENSION(40) :: D=0.0 ! implied/implicit loops for initial values
The actual allocation of memory for dynamically allocatable arrays must be done within the executable code after the declarations:
... NSTEP=1000 ... ALLOCATE(X(NSTEP),Y(NSTEP)) ! allocate array space as executable ...
I=?? ! get integer subscript from somewhere X(I)=5.0 ! reference array element using integer Y(I)=X(I)+10.0
READ *, R1,R2 ! use the default input stream PRINT *, "The sum is ",SUM ! use the default output stream
OPEN(UNIT=1,FILE="flight.dat",STATUS="OLD") ! open an existing file READ(UNIT=1,FMT=*) T,H ! read from a unit number ... OPEN(UNIT=2,FILE="output.dat",STATUS="NEW") ! open a new file WRITE(UNIT=2,FMT=*) T,H ! write to a unit number ... CLOSE(UNIT=1) ! close a unit number
PRINT 10, II,RR ! format specified by a labelled statement
10 FORMAT(5X,I3,1X,F10.3)
WRITE(2,"(5X,I3,1X,F10.3)") II,RR ! format specified in-line
The following field descriptors are available:
nX n spaces
Fn.m REAL number, fixed point, m digits after point
En.m REAL number with exponent, m digits after the point
In INTEGER
An CHARACTER string
Arithmetic operators:
( ) parentheses to dictate the order of operations ** exponentiation ( raising to a power) * or / multiplication or division + or - addition or subtraction
INTEGER :: J REAL : X,Y,A,B ... X=A+B ! set variable X equal to A+B J=J-1 ! decrement variable J by 1 Y=X**2 ! set variable Y to X squared
INTEGER :: I,L,II REAL : A,R,R1,R2 ... I=R+35.3 ! The REAL result is truncated before assignment to INTEGER A=I/L ! The result is truncated to an INTEGER and assigned to REAL R2=R1+II ! The result of mixed arithmetic is REAL
The following relational operators are available:
== or .EQ. equals /= or .NE. not equals > or .GT. greater than < or .LT. less than >= or .GE. greater or equal <= or .LE. less or equal
The following logical operators are available:
.NOT. invert logical value .AND. logical and .OR. logical or .EQV. logical equivalence .NEQV. logical non-equivalence
Decision structures:
IF(A<3.0) THEN ! start of an IF structure ... ELSE IF(A>10.0) THEN ! continue IF structure ... ELSE ! final clause of IF structure ... END IF ! end of IF structure
SELECT CASE(expression_select) ! start of a CASE structure
CASE(expression1) ! expression CHARACTER or INTEGER
! do something
CASE(expression2)
! do something else
CASE DEFAULT
! some default
END SELECT ! end of CASE structure
Definite loops:
DO J=1,10 ! start a definite loop with counter J
! within loop J will be 1,2,....9,10
END DO
DO K=1,13,2 ! start definite loop with increment 2
! within loop K will be 1,3,5....9,11,13
END DO
Indefinite loops:
NMAX=100 ! set upper limit
NMIN=-100 ! set lower limit
DO ! start indefinite loop
J=?? ! get J from somewhere, somehow!
IF((J>NMAX).OR.(J<NMIN)) THEN ! test for exit
EXIT ! jump to after END DO
ELSE
! keep going
ENDIF
END DO
...
Implied loops:
REAL , DIMENSION(10) :: A ... READ(UNIT=1,FMT=*) (A(K),K=1,10) ! implied loop on input
implicit loops for arrays:
REAL, DIMENSION(100) :: A,B,C ... A=1.4 ! set all elements of array A to a scalar B=A+4.0 ! add 4.0 to all elements of A C=A*B ! multiply elements of A by elements of B D=SIN(C) ! take the SIN of each element of C
ABS absolute value SQRT square root SIN sine (argument radians) COS cosine (argument radians) TAN tangent (argument radians) ASIN arcsine (returns radians) ACOS arccosine (returns radians) ATAN arctangent (returns radians) ATAN2(X,Y) arctangent Y/X (returns radians) ALOG10 logarithm base 10 EXP raise e to power INT truncate REAL to INTEGER NINT fix REAL to nearest INTEGER REAL float INTEGER to a REAL MOD remainder of A1/A2
A single PROGRAM module provides the entry point from the operating system:
PROGRAM progname ... END PROGRAM progname
SUBROUTINEs can be called from within PROGRAM or other SUBROUTINEs or FUNCTIONs:
SUBROUTINE subname(dummy argument list) ... END SUBROUTINE subname
FUNCTIONs can be used within PROGRAM, SUBROUTINEs, or other FUNCTIONs:
FUNCTION funname(dummy argument list) ... END FUNCTION subname
A MODULE is used for the declaration of variables and parameters that are to be shared (global) between PROGRAM, SUBROUTINEs and FUNCTIONs:
MODULE modname
declarations...
END MODULE modname
SUBROUTINE whatever
USE modname ! use declarations in a MODULE
...
Calling subroutines and using functions:
SUBROUTINE DOIT(A,B)
IMPLICIT NONE
REAL :: A,B ! declare dummy arguments
internal declarations
executable instructions
END SUBROUTINE DOIT
FUNCTION ZAP(C)
IMPLICIT NONE
REAL :: ZAP,C !declare function and dummy arguments
internal declarations
executable instructions
ZAP=expression !assign value to function before return
END FUNCTION ZAP
PROGRAM TRY
REAL :: X,Y,Z
X=4.0
CALL DOIT(X,Y)
Z=ZAP(X)
PRINT *, X,Y,Z
END PROGRAM TRY
A safe way of passing arrays as arguments:
SUBROUTINE THINGY(N,A,B)
IMPLICIT NONE
INTEGER :: N
REAL, DIMENSION(N) :: A
REAL :: B
...
END SUBROUTINE THINGY
PROGRAM TRYIT
INTEGER, PARAMETER :: NEL=10
REAL :: X(NEL),Y(NEL)
...
CALL THINGY(NEL,X,Y)
...
END PROGRAM TRYIT