Document Sample

Chapter 1 Programming in a high-level language – Fortran 90/95 Programming Class 1 (Week 1) Why learn a high-level language? The cutting edge in scientiﬁc computing is eﬃcient implementation of ﬂoat- ing point arithmetic, i.e. arithmetic involving numbers expressed in decimal expansions. Although integrated systems such as Matlab are highly convenient you pay a penalty for this in terms of eﬃciency and memory requirements. Matlab on its own is not suitable for very large problems (or even medium-sized ones without careful coding). High-level compiled languages (Fortran, C, C++, Pascal) are more eﬃcient. Which high-level language to learn? Fortran77. “Old fashioned”, for example it does not allow abstract data types or recursion. However, since a great deal of well-engineered scientiﬁc software is in Fortran77, it remains important today. The simplicity of its syntax al- lows compilers to highly optimise the resulting machine code making full use of the cache, the pipeline and other hardware features in modern computers. Fortran90/95. modern upgrade of Fortran77 which supports many important fea- tures missing from Fortran77. It is backward compatible with Fortran77. (Fortran2003 standard coming!) C and C++. Many scientiﬁc programmers also use languages like C (especially) 1 or C++. In terms of eﬃciency comparable to Fortran. Pascal. “Old fashioned” like Fortran77. Java. Not as eﬃcient for scientiﬁc applications. There is a growing interest in object-oriented programming languages (like C++ or Java) for scientiﬁc applications as opposed to procedural languages (like Fortran). Modern C++ compilers allow for similar eﬃciency as Fortran77, while allowing for data hiding, code abstraction and other modern developments in computer science. (Fortran2003 standard includes more object-orientation!) Many workers in scientiﬁc computing will need to work with multi-language codes (and indeed it is nowadays quite easy to couple Fortran and C or C++ code). However, for us it is more important to learn the principles of software design, then to learn a particular language. In this course we will use Fortran95 and we will spend the ﬁrst 2 weeks learning the basics of this language. Example 1 Assemble the n × n Hilbert matrix Hi,j = 1/(i + j − 1) , i, j = 1, . . . n . Programs given: • Fortran95 program in week1/hilbert.f90 • Three Matlab programs in week1: hilbert.m, hilbert1.m and hilbert2.m Lab Exercise 1. Connect to the computers in the MSc room and edit these ﬁles. 1. Start → Programs → Unix Applications → X-Windows Manager 2. On grey area right click, select any BUCS machine. Logon. 3. Type ssh -l<userid> mapc-desknode-0**.maths 4. Copy my ﬁles over to your ﬁle space: cp ˜ jvl20/ma50177/week1/hilbert* . 5. Type emacs& to start the Emacs-editor. Use the file button to open hilbert.f90. 2 The Fortran95 source code in week1/hilbert.f90: !------------------------------------------------------------------ program hilbert ! optional program label ! f95 program to assemble the Hilbert matrix ! We begin by declaring our variables real, allocatable, dimension(:,:) :: H ! declares a real matrix ! H with two dimensions ! (i.e. rank 2). The extent ! of each dimension is declared ! in the allocate command below real :: t1,t2 ! real numbers used for timing ! this program integer :: i,j,n ! integers, used for counting print*,’Type n’ ! message for the user read*, n ! extent of the dimensions of H allocate(H(n,n)) ! allocate storage call cpu_time(t1) ! look at the clock do i = 1,n ! nested loop for assembling H do j = 1,n H(i,j) = 1.0/(i+j-1.0) ! mixed mode expression for H(i,j) end do end do call cpu_time(t2) print*,’n=’,n,’cputime taken = ’,t2-t1 ! prints out timing results print*, H(1,2) ! prints out (1,2)th element of H end program hilbert ! mandatory end statement ------------------------------------------------------------------- 3 To compile this program type ifort hilbert.f90 to the $ prompt. You will get an executable a.out . Type a.out to run this executable. Notes on hilbert.f90: 1. “Dynamic” allocation of array storage for H. 2. All variables must be declared! 3. Study the syntax of the do loop. 4. Real and integer data types appear – diﬀerent in Fortran95 ! 5. “Mixed mode” arithmetic in the formula for H(i,j) ([MeRe p.36], [Sm p.16]). 6. Use of intrinsic function cpu_time. 7. Only H(1,2) printed to screen, as example. Lab Exercise 2. What is the eﬀect of the following changes to hilbert.f90? 1. Change H to h in the last print statement. 2. Comment out the line: real,allocatable... (use ! to comment out a line). 3. Comment out the line: integer :: i,j,n 4. Change the line integer :: i,j,n to integer :: i,j,n,k Make each of the changes, recompile and rerun to see the result. What do you learn from this? Now, repeat these experiments using the compile command: ifort -u . Try some other changes. To illustrate the eﬃciency of a compiled language, compare the program hilbert.f90 with two Matlab programs which carry out the same task: 4 %------------------------------------------------------------------ % program hilbert.m : script to assemble the n x n Hilbert matrix % using a double loop n = input(’type n: ’) H = zeros(n,n); % Initialise H as an n x n matrix tic % look at clock for j = 1:n % Use a double do loop for i = 1:n H(i,j) = 1/(i+j-1); % doesn’t matter whether we % write 1 as 1.0 or not end % end of j loop end % end of i loop toc H(1,2) % end of program hilbert.m %------------------------------------------------------------------ %------------------------------------------------------------------ % program hilbert2.m : script to assemble the n x n Hilbert matrix % avoiding loops by using array manipulation facilities of Matlab n = input(’type n: ’) H = zeros(n,n); % Initialise H tic H = ones(n,n)./((1:n)’*ones(1,n) + ones(n,1)*(1:n) - ones(n,n)); toc H(1,2) % end of program hilbert2.m %------------------------------------------------------------------ 5 Eﬃciency: CPUtime in seconds (2.1GHz Athlon, 3.7 GByte RAM) Matlab 6 (Release 13) Matlab 7.3 (Release 2006b) n hilbert.m hilbert2.m hilbert.m hilbert1.m hilbert2.m hilbert.f90 1000 5.4 0.29 0.038 0.067 0.14 0.018 2000 21.6 1.13 0.15 0.61 0.53 0.075 4000 88.6 4.5 0.60 4.2 2.1 0.30 8000 **** **** 2.4 17.1 8.4 1.2 16000 **** **** **** **** **** 4.9 (**** = out of memory (3.7 GByte RAM)) In Matlab 6 (and in any older version) hilbert.m is almost 300 times slower than the compiled code hilbert.f90 and hilbert2.m is still 15 times slower. In Matlab 7 the times are closer, but the compiled code is still more than twice as fast. The high cost of hilbert.m (in Matlab 6) is due to the fact that Matlab is an interpreted language: each line of code is translated to machine instructions each time it is reached: here the same line is retranslated n2 times. We can see that in the later versions of Matlab (since Release 14, September 2004) the interpreter is able to identify the loop and to avoid retranslating the same line n2 times. However, this will not always be possible when a for-loop becomes more complicated. The faster code hilbert2.m relies on Matlab built-in matrix manipulation functions to avoid loops. This built-in functions are themselves Fortran routines from fast libraries like BLAS and LAPACK which we will come across later in this course. A further problem with interpreted languages like Matlab is memory. Even on my computer, which has 3.7 GByte memory, I can only go with n up to 4000 or 8000 for Matlab 6 and Matlab 7.3 respectively, whereas in hilbert.f90 I can go up to 20000! The bottom line is that Matlab has problems to compete for large problems. This will become even more apparent for more complicated problems. 6 Lab Exercise 3. Complete the two tables on any of the linux machines in the MSc room (3.0Ghz Intel, 1 GByte memory, Matlab 7 (Release 14)). n hilbert.m hilbert1.m hilbert2.m hilbert.f90 1000 2000 4000 8000 16000 7 Programming Class 2 (Week 1) IMPORTANT: You should read this already before your Thursday class. The exercises will be done in the class. Please come prepared to ask questions. This class will take place in the MSc room. The purpose is: • to learn some Emacs skills • to learn about elementary arithmetic, data types, and intrinsic functions in Fortran95 • to learn about control of ﬂow in Fortran95 programs • to translate an elementary MATLAB program into Fortran95 By the end of Friday’s class you should have done all the exercises in this section. The Emacs editor Emacs is the recommended editor for developing Fortran95 code. To start Emacs type emacs& to the $ prompt. You can use the menu bar at the top to load, insert and save ﬁles. You can use Emacs as a ﬁle manager. Under the ‘ﬁle” button, select “Open direc- tory”. You can then navigate around your ﬁles and if you like you can use Emacs to create and delete directories and ﬁles. To get pretty colours for your Fortran code, go to the Options menu and select Syntax Highlighting. Then click on Save Options. Compiling and running programs Exercise 1. A simple program for numerical integration. Suppose f is a given function on an interval [α, β]. Then Simpson’s rule for com- β puting an approximation to the exact integral α f (x)dx is: h Iα,β (f ) = {f (α) + 4f (γ) + f (β)} , 6 where h = (β − α) and γ = (α + β)/2 . A program to implement this in a very simple way is at 8 ~jvl20/ma50177/week1/s_simpson.f90 Copy this program to your ﬁlespace. Compile and run it. What does the result tell you? Replace f (x) = x3 by f (x) = x5 . What do you ﬁnd? Intrinsic functions and data types Exercise 2. Intrinsic functions. Look up the list of intrinsic functions in the Fortran95 book of Smith [Sm] or Met- calf and Reid [MR]. Find the entry for atan. Using the fact that tan(π/4) = 1, write 1 a line of code to compute π and then compute an approximation to 0 π cos(πx)dx using Simpson’s rule. What is the error in the result? Exercise 3. Data types. Fortran95 allows reals of both single and double precision. You can specify the precision by specifying a kind parameter: • real (kind=4) ... a single precision real • real (kind=8) ... a double precision real Note that if you do not specify the kind parameter, single precision is used by default. Also, if a program requires the real constant 10.0 in double precision, it should be typed as 10.0_8. Otherwise it will only be of single precision and might ﬂaw any double precision calculations. Edit your s_simpson.f90 so that the computations are done in double precision. Compile and run again. How can you be sure that the double precision is working? Hint: To avoid having to edit the declaration of every variable when you want to change the precision of a program, you can simply declare the kind as a parameter and then you only need to change it once. So in the declaration of variables you add something like: integer, parameter :: rk=1 Then when declaring a real variable integral you say real (kind=rk) :: integral 9 Also, if a program requires the real constant 10.0, it can be typed as 10.0_rk to obtain the required precision. Then when editing you only need to change the value in the parameter declaration once. Warning: Fortran95 is a weakly typed language. Variables that start with i,j,k,l,m or n are understood implictly to be integers unless declared otherwise. This is a rather dangerous Fortran77 legacy. To override it you can put implicit none at the top of each program unit. Alternatively, compile with the -u option. Recom- mendation: always compile with the -u option. Writing your ﬁrst program Exercise 4. Control of Flow. As you will remember in Matlab , the most common constructions for controlling the ﬂow in a program are the (i) for loop, (ii) while loop, (iii) if statement. In Fortran95 the counterpart of (i) is a do loop. You have seen an example of this already in hilbert.f90 . For (iii) there is an if statement which is similar to Matlab but with slightly diﬀerent syntax. For (ii), there is no explicit while loop. However, the required action can be achieved with a do loop combined with the if and exit statements. For example, to loop over all values of i from 1 to n, stopping when i=j we can write do i = 1,n . . if (i==j) exit . . end do 10 The exit command causes exit from the do loop and the program continues with the command after the end do statement. Look in the Fortran95 books for more detail. What to do: Copy over the program ~jvl20/ma50177/week1/macheps.m. This is a Matlab program to compute the machine epsilon for a given machine. You should have seen this before in MA50174. The machine epsilon is the smallest positive number ǫ such that the machine representation of 1 + ǫ diﬀers from the machine representation of 1. Translate this program from Matlab to Fortran95. Compile and run your Fortran95 code and compare its output to the Matlab code. IMPORTANT: Please try to complete these exercises before Friday’s class. 11 Programming Class 3 (Week 1) SUMMARY: Main Program and Subprograms Subprograms (Procedures): function and subroutine The argument list, intent of an argument Interface blocks The Makeﬁle Examples: various quadrature procedures Example: Numerical Integration β Simpson’s rule for approximating α f (x)dx is h Iα,β (f ) = {f (α) + 4f (γ) + f (β)} 6 where h = (β − α) and γ = (α + β)/2 . Task 1. Open the ﬁle s_simpson.f90 in emacs (see also below). !-------------------------------------------------------------------- program s_simpson ! simple program to approximate the integral of a function f over the ! interval [alpha, beta] using Simpson’s rule. Here f(x) = x**3 real :: alpha, beta, gamma real :: f_alpha, f_beta, f_gamma, h real :: integral print*, "Type alpha,beta" read*, alpha,beta gamma = (alpha + beta)/2 12 h = (beta-alpha) f_alpha = alpha**3 f_beta = beta**3 f_gamma = gamma**3 integral = (h/6)*(f_alpha + 4*f_gamma + f_beta) print*, integral end program s_simpson !-------------------------------------------------------------------- It works but is not very convenient/ﬂexible: 1. To change the function f we have to edit 3 lines and recompile. 2. We often will want to apply the rule over many small subintervals of a larger interval. This is the “composite Simpson’s rule”. To do this we would have to rerun the current program many times, type in new data each time, and add up the answers afterwards. To ﬁx Problem 1 we put the deﬁnition of f in a function subprogram: function func (x) result (func result) real, intent(in) :: x real :: func result func result = x**3 end function func This is stored as func.f90 (in week1/simpson). The intent statement is optional here. It is good programming practice to put it in. x is a dummy argument in func. If a dummy argument is speciﬁed as intent(in) then it cannot be changed inside the function. Now in s_simpson we can write f_alpha = func(alpha) Then alpha is the actual argument, which is passed to func. 13 To address Problem 2, it would be better to rewrite the whole Simpson procedure as a subprogram, which takes as arguments the function func and the end points alpha, beta and produces the approximate integral. This could then be called many times from a main program, with diﬀerent values of the arguments, if required. The best way to do this is using a subroutine. Unlike a function, a subroutine often has both input and output arguments, or even arguments which have intent(inout). It has no “result” as such. All data transfer with the calling program is via the arguments. Also, unlike a function, a subroutine is invoked with the call command. !-------------------------------------------------------------------- subroutine simpson(integral,func,alpha,beta) ! simple program to approximate the integral of a function func ! over the interval [alpha, beta] using Simpson’s rule real, intent(out) :: integral ! The result of Simpson’s rule real, intent(in) :: alpha, beta ! The inputs real :: gamma, h ! variables needed internally ! Interface block to ensure consistency and detect errors in passing ! arguments to subprograms essentially to let the subroutine know ! what type of thing func is interface function func(x) result(func_result) real, intent(in) :: x real :: func_result end function func end interface gamma = (alpha + beta)/2 h = (beta-alpha) integral = (h/6)*(func(alpha) + 4*func(gamma) + func(beta)) end subroutine simpson !-------------------------------------------------------------------- This subroutine is called from a main program, in our case integration.f90 14 All these ﬁles are located in week1/simpson. Task 2. Copy the entire directory week1/simpson by typing: cp -R ~jvl20/ma50177/week1/simpson . (‘‘-R’’ for recursive) Task 3. Compile with the command ifort -u func.f90 simpson.f90 integration.f90 Code development: The Makeﬁle Since we now have three program units which can be developed separately it is best to control their development using a Makeﬁle. This contains all the commands which are needed to compile the program units and then make an executable. When you have a Makeﬁle, type the command make to run it. (Only program units which have been edited/changed since the last time you ran make are recompiled.) A suitable Makeﬁle to compile and link all the components of our code in the ﬁles func.f90 simpson.f90 integration.f90 is in week1/simpson/Makefile: #-------------------------------------------------------------------- # This is a comment line in the Makefile OPT = -u # all ifort calls use the -u option # (suppresses implicit typing) all: integration # forces all the commands needed to create # ‘‘integration’’ to be executed, # unless there is nothing to do func.o: func.f90 # compiles func.f90. The -c option ifort $(OPT) -c func.f90 # causes the production of an object # file only which will be linked # with other object files below # then does the same for simpson.f90 and integration.f90 simpson.o: simpson.f90 ifort $(OPT) -c simpson.f90 15 integration.o: integration.f90 ifort $(OPT) -c integration.f90 # now link the three object files to get an executable # -- in this case called integration integration: integration.o simpson.o func.o ifort -o integration $(OPT) integration.o simpson.o func.o clean: rm -rf *.o core # creates a command to remove old .o and # core files # type ‘‘make clean’’ to run this #-------------------------------------------------------------------- Task 4. Go to the directory simpson in your ﬁlespace and type make clean and then make Note. Using ifort -o testx test.f90 creates and executable called testx (rather than a.out). Passing parameters Let us return to the example of approximating 1 π cos(πx)dx. 0 To implement Simpson’s rule we can edit the ﬁle func.f90 to compute π cos(πx). However we need to evaluate π. It would be wasteful to have the function func evaluate π every time it is called. It is better to evaluate π once (in the main program) and then pass it down, ﬁrst to the simpson subroutine and then to func. A suitable implementation of this is in week1/simpson2 Task 5. Copy over the whole directory with the command: cp -r ~jvl20/ma50177/week1/simpson2 . Have a look at the ﬁles simpson2.f90, func2.f90, integration2.f90 16 Lab Exercise 1 Create a suitable Makeﬁle in your directory simpson2. Check that it works. Lab Exercise 2 (Composite Simpson’s rule) The 3 point Simpson’s rule is not accurate enough for most integation tasks. To im- prove accuracy, the interval of integration [α, β] may be divided into N subintervals [xi−1 , xi ], where xi = α + i(β − α)/n, i = 1, · · · , N . Write a subroutine csimpson with input arguments func, α, β and N which implements this, by calling simpson on each subinterval. Test it by integrating func(x) = ex , over [0, 1], and by checking that the approx- imation converges to the exact integral as N → ∞. 17

DOCUMENT INFO

Shared By:

Categories:

Tags:
fortran 90 95, fortran 90, stephen chapman, how to, customer reviews, fortran 95, sun grid engine, batch jobs, programming languages, cluster nodes, parallel job, stephen j. chapman, intel cpus, free ebook, programming language

Stats:

views: | 32 |

posted: | 2/1/2010 |

language: | English |

pages: | 17 |

OTHER DOCS BY hcw25539

How are you planning on using Docstoc?
BUSINESS
PERSONAL

By registering with docstoc.com you agree to our
privacy policy and
terms of service, and to receive content and offer notifications.

Docstoc is the premier online destination to start and grow small businesses. It hosts the best quality and widest selection of professional documents (over 20 million) and resources including expert videos, articles and productivity tools to make every small business better.

Search or Browse for any specific document or resource you need for your business. Or explore our curated resources for Starting a Business, Growing a Business or for Professional Development.

Feel free to Contact Us with any questions you might have.