Fortran

Markus Meuwly Michael Devereux

Stefan Willitsch

Department of Chemistry University of Basel

Einführung in die Programmierung mit Fortran

Fortran

Ziel / Motivation

Sie sollten in der Lage sein: • Einfache Programme zur Lösung naturwissenschaftlicher Probleme zu

schreiben, zu kompilieren und zu benutzen.

Weshalb Programmieren? • Geschwindigkeit • Genauigkeit • Flexibilität

Fortran

übersicht

• Einführung • Auswertung von Funktionen (function, subroutine) • Numerische Differentiation • Numerische Integration • Matrixoperationen

Fortran

übersicht

Die Programmiersprache Fortran ist eine sog. “höhere Programmiersprache” im Gegensatz zu Assembler. Es gibt von Fortran verschiedene “Dialekte” Fortran I, II, IV Fortran 66 Fortran 77 Fortran 90, 95, 2000, 2003, 2008, 2010 Wir werden uns mit Fortran 77 beschäftigen. Neben Fortran gibt es noch C, C++, welche im high performance computing verwendet werden. Allen höheren Programmiersprachen gemein ist die Tatsache, dass ein Computerprogramm in der entsprechenden Sprache geschrieben wird und dann durch einen Compiler in ein ausführbares Programm übersetzt wird.

Fortran

Allgemeine Programmstruktur

program test !Name des Programms

Implicit none !Ausschalten der autom. Type conversion

Integer i !Deklaration einer integer-Zahl

Real*8 r !Deklaration einer reellen Zahl

Complex cn !Deklaration einer komplexen Zahl

Character c1 !Deklaration eines alphanumerischen String

Etc…

Input !Eingabe/Einlesen von Daten

Instructions !Anweisungen

Output !Ausgabe der Resultate

end !Ende des Programms

Fortran

Allgemeine Programmstruktur

Einige Merkmale von Fortran: Gute Programmierpraxis ist es, die automatische Typenkonversion auszuschalten. So vermeidet man eine Vielzahl von Fehlern. Also: implicit none In Fortran 77 beginnen die Instruktionen in der 7. Spalte des Editors – manche Editoren erkennen dies selbständig. In späteren Fortran Versionen fällt dieser Umstand weg. Die erste Spalte kann für Kommentare verwendet werden: C This is a comment and will not be interpreted by the compiler

Die maximale Zeilenlänge sind 80 Zeichen; ist ein Kommando länger, muss es umgebrochen und fortgesetzt werden; das Fortsetzungszeichen “&” steht in der 6. Spalte: x = tanh(sin(x))-dble(i)*erf(x)*dexp(-sin(x))*

& x**2-4*x

Fortran unterscheidet nicht zwischen Gross- und Kleinschreibung

Fortran

Eingabe/Ausgabe Anweisungen

Format: Read(nunit,format)

Write(nunit,format)

Open(nunit,name)

Close(nunit)

nunit ist eine Zahl, welche es erlaubt, ein bestimmtes File anzusprechen (s. Open) Format ist eine Formatierungsanweisung, welche es z. B. erlaubt, eine gewisse Anzahl Stellen bei der Ausgabe von Resultaten anzuzeigen. Beispiele: Open(unit=3,file=‘name.dat’)

Read(3,*) zahl

Read(3,100) zahl

100 format(f12.5)

C F12.5 steht fuer floating point (reelle Zahl)

C mit insgesamt 12 Stellen und 5 Stellen nach dem Komma

write(22,101) zahl

101 format(f10.4)

Fortran

Kompilieren und Ausführen von Programmen

Zweck ist es, aus dem Quellcode (“lesbar”) ein ausführbares Programm (“nicht-lesbar”) zu machen. Das “executable” run wird aus dem Quellcode run.f folgendermassen erzeugt:

>gfortran –o run run.f

Dabei steht “>” für den Cursor unter UNIX. Von der Kommandozeile aus wird das Programm folgendermassen ausgeführt (”laufen gelassen”): >./run

Benötigt das Programm einen “input” – also sollen Daten extern eingelesen werden – so geschieht dies mittels einer “Pipe”. Dabei wird der input input.dat mit “<“ hinein-gepiped und der output mit “>” in eine output-datei umgeleitet (statt auf Bildschirm geschrieben). >./run < input.dat > output.dat

Aufgabe: Programm zur Berechnung der Flaeche eines Kreises mit Eingabe des Radius.

Fortran

Einfache Funktionen und Prozeduren Sobald Programme etwas umfangreicher werden, empfiehlt es sich, gewisse Teile in Funktionen und Unterprogramme auszulagern. Die erlaubt es auch, wiederkehrende Aufgaben effizient zu erledigen. Bei einer function wird am Ende dem Funktionsnamen ein entsprechender Wert zugewiesen Beispiele: Suchen von Nullstellen mittels Newton-Raphson Diagonalisieren von Matrizen

Die allgemeine Struktur einer function function f(…) ! Name der function muss mit Wert am Ende übereinstimmen

Deklarationen ! Wieder mit implicit none etc.

Anweisungen !argumente werden verarbeitet

f=…

…

return

end

Fortran

Einfache Funktionen und Prozeduren Im Gegensatz zu einer function wird dem Namen einer subroutine kein Wert zugewiesen. Vielmehr werden zur Ein- und Ausgabe die Argumente der subroutine verwendet. Die allgemeine Struktur einer subroutine

subroutine name(…) !name sollte die Aufgabe/Methode widerspiegeln

Deklarationen

Anweisungen

…

return

end

Beispiele sind Diagonalisieren von Matrizen, Ordnen einer Liste Aufgaben Programm zur Auswertung einer algebraischen Funktion Programm zur Nullstellenberechnung

Fortran

Einfache Funktionen und Prozeduren Beim Auswerten von Funktionen ist auf numerische Instabilitäten zu achten. Zum Beispiel ist die Auswertung des Ausdrucks für grosse x unstabil, da für grosse x die Differenz zwischen x und x-1 sehr klein wird; dementsprechend ginge das Resultat gegen 0, was jedoch nicht sein kann, was eine Taylorentwicklung zeigt: Die numerisch stabile Variante ist: Diese Probleme werden unter dem Begriff “numerische Auslöschung” oder “cancellation” zusammengefasst. Weitere Beispiele, die durch verschiedene mathematische Umformungen in stabile Varianten übergeführt werden können sind: Für kleine x; entwickle sin(x) in Taylorreihe Für kleines e; forme mit Trigonometrie um Für grosses x; erweitern – siehe oben

xx −+1

( ) ( )

( )( )( ) xxxx

xxxxxx

xx

xxxx

xx

++=

++++

−+=−+

+−=++−=

+−−−+−

+−+

≈−+

11

1111

41

85

81

41

21

1811

211

81

211

1

22

....

....

( )

111+

−

−+

xx

xxx

x

sinsin

sin

ε

Fortran

Einfache Funktionen und Prozeduren: Finden von Nullstellen Finden von Nullstellen mittels • Bisection • Regula Falsi • Newton Raphson

Bei den ersten beiden Methoden nutzt man aus, dass das Produkt der Funktionswerte das Vorzeichen ändern muss, wenn dazwischen eine Nullstelle ist. Bisection: Dann wird ein dritter Punkt zwischen x1 und x2 definiert und die Vorzeichenänderung verfolgt. Auf diese Art kann das Intervall immer enger eingeschränkt. Regula Falsi: Hier wird eine lineare Interpolation gemacht, für welche ein neuer, verbesserter unterer Wert analytisch bestimmt werden kann: Die beiden ersten Methoden unterscheiden sich dadurch, wie rasch sie konvergieren. Newton Raphson

( )( ) ( )( ) ( ) 0

00

21

21

21

<<>

xfxfxfxf

xx;

,

( )( ) ( )

( ) ( )( ) ( )12

12211

21

21

00

xfxfxfxxfxx

xfxfxx

−−

=

><

~

;,

( )( )k

kkk xf

xfxx′

−=+1

Fortran

Einfache Funktionen und Prozeduren FUNCTION rtnewt(x0,xacc) ! Newton-Raphson Funktion

implicit none

real*8 x0,xacc ! Deklaration der externen Variablen

integer JMAX,j ! Deklaration der internen Variablen

real*8 rtnewt,df,dx,f

external funcd ! Deklaration der Funktions-Subroutine

parameter (JMAX=50) ! Def. Parameters; Anzahl Iterationen

rtnewt=x0 ! Start-Wert

do j=1,JMAX ! DO-LOOP (root-finding)

print*,'cycle ',j,'r_trial=',rtnewt

call funcd(rtnewt,f,df) ! Call Morse-Subroutine 'funcd'

dx=f/df

rtnewt=rtnewt-dx

if(abs(dx).lt.xacc) return ! Convergence

enddo ! Ende des Loops

print*,'rtnewt exceeded maximium iterations (',JMAX,')'

END ! Ende der Funktion

Fortran

Einfache Funktionen und Prozeduren program main

implicit none

Deklarationen

Eingabe

Anweisungen ! Incl. calls to subroutines etc.

Ausgabe

End

C--------------------

Subroutine sb1(arguments)

implicit none

Deklarationen

Eingabe !Oft unnötig, da in “arguments” übergeben

Anweisungen

Ausgabe !Oft nur als test

End

C---------------------

Function fct1(arguments)

implicit none

Deklarationen

Eingabe !Oft unnötig, da in “arguments” übergeben

Anweisungen

Ausgabe !Oft nur als test

C--------------------------

Weitere subroutines und functions

Fortran

Numerische Ableitung und Integration

Sind Funktionen nur auf einem diskreten Gitter bekannt (beispielsweise durch Messen einer Grösse als Funktion der Temperatur welche nicht kontinuierlich reguliert werden kann), so stellt sich das Problem, Ableitungen numerisch zu berechnen. Die dazu üblichen Formeln heissen “finite differences”. Wenn man die einfache Definition einer ersten Ableitung unkritisch anwendet, können die Ableitungen sehr ungenau werden. 2 Fehlerquellen führen dazu Abschneidefehler Rundungsfehler Ein stabiler Algorithmus wird in den Uebungen vorgestellt. Aufgabe: Programm zur numerischen Ableitung einer algebraischen Funktion Newton-Raphson mit numerischer statt analytischer Ableitung

Fortran

Numerische Ableitung und Integration

( ) ( ) ( )h

xfhxfxf −+≈′

Probleme: Wenn man blind einen Wert für h einsetzt, ergeben sich numerische Instabilitäten, da sich x und x+h nicht beliebig genau im Rechner darstellen lassen. Es ist sicherer, bei gegebenem x, die Schrittweite h durch folgende Konstruktion zu bestimmen: tmp = x+h h=tmp-x Durch Verwenden einer symmetrisierten Formel für die Ableitung erhält man ebenfalls eine Stabilisierung:

( ) ( ) ( )h

hxfhxfxf2

−−+≈′

Fortran

Numerische Ableitung und Integration

Numerische Integration ist eine wichtige Disziplin in der numerischen Behandlung von naturwisenschaftlichen Problemen. Man unterscheidet dabei zwei Kategorien: Die Integrationspunkte sind regelmässig verteilt Numerische Quadraturen – die Punkte sind nicht gleichmässig verteilt. In allen Fällen wird ein Integral durch eine Summe approximiert. Auf einem äquidistanten Gitter xi=x0+ih ergeben sich die Trapezregel und die Simpson Regel, welche Polynome ersten, zweiten und dritten Grades exakt integrieren. Für andere Funktionen ergeben diese Formeln einen Näherungswert für das Integral.

( )

( )∫

∫

++≈

+≈

3

1

2

1

3312

341

31

2211

21

x

x

x

x

fffhdxxf

ffhdxxf

Fortran

Numerische Ableitung und Integration

Für die Romberg Integration wird der Fehlerterm in der Simpson-Regel explizit behandelt, und zwar durch eine Polynominterpolation. Aufgabe: Programm, welches nacheinander das Integral einer Funktion ueber ein endliches Intervall mit der Trapezregel, der Simpson Regel und der Romberg Integration berechnet. Hinweise zur Lösung der Aufgaben: Sie erhalten mehrere subroutines trapzd, qtrap, qsimp, qromb, polint

Diese können Sie unverändert verwenden. Beim Kompilieren müssen diese Programme mit an das Hauptprogramm gelinkt werden: gfortran –o prog main.f trapzd.f etc….

Dabei ist main.f Ihr Hauptprogramm.

Fortran

Numerische Ableitung und Integration program main

implicit none

Real*8 func,a,b,s

a=xx

b=xx

call qsimp(func,a,b,s)

write(6,*) ‘Wert des Integrals ‘,s

end

C--------------------

function func(x)

implicit none

real*8 func,x

func=xxxx

return

end

Fortran

Matrixoperationen

Matrizenrechnung ist eine weitere wichtige Anwendung. Dabei stehen folgende Operationen im Vordergrund: Matrixmultiplikation Lösen linearer Gleichungssysteme Diagonalisierung (Eigenwerte, Eigenvektoren) Aufgaben Programm zur Matrixmultiplikation Lösen eines linearen Gleichungssystems (Gauss Elimination)

Fortran

Matrixoperationen Eine Matrix ist ein rechteckiges Schema welches reelle oder imaginäre Zahlen enthält. In Fortran wird eine Matrix folgendermassen deklariert: Real*8 a(2,2)

Dieses Statement deklariert eine 2x2 Matrix. Dieses Schema entspricht der folgenden Anordnung: Computertechnisch heissen Matrizen allgemeiner “arrays”. Sie können auch mehr als 2 Indizes aufweisen: Real*8 a(2,2,4,2)

Ausdrücke, welche in Matrixform geschrieben werden, sind Ihnen beim Lösen von Gleichungssystemen begegnet:

Beachten Sie, dass die Indizierung von arrays von “1” losgeht. Das erste Element in einem array ist also a11 und nicht a00. Sollte es notwendig sein, kann man die Indizes auch von “0” weg laufen lassen: Real*8 a(0:2,0:2)

( )2141

1235

42135

=−=

=

−

=+−=+

yxyx

yxyx

,

2221

1211

aaaa

Fortran

Matrixoperationen Nachdem die Matrix deklariert wurde, muss sie mit Werten gefüllt werden. Dies geschieht auf die übliche Weise einer Zuordnung. Ein kleines Beispielprogramm: Program matrix

Implicit none

Real*8 a(2,2)

A(1,1) = 1.d0

A(1,2) = 0.d0

A(2,1) = 0.d0

A(2,2) = 1.d0

end

Dieses Fragment definiert eine 2x2 Einheitsmatrix. I.a. sind für diesen speziellen Fall “do-loops” vorgesehen:

Fortran

Matrixoperationen Program matrix

Implicit none

Real*8 a(2,2)

Integer I,j

Do I = 1,2

Do j = 1,2

A(I,j) = 0.d0

If (I .eq. j) a(I,j) = 1.d0

Enddo

Enddo

end

In der Uebungsaufgabe sollen Sie ein Gleichungssystem lösen. Dazu werden Ihnen wieder einige subroutines zur Verfügung gestellt. Die Aufgabe besteht demzufolge darin, die Matrix A sowie den Lösungsvektor b mit entsprechenden Werten zu füllen und diese an die subroutine zu übergeben, welche das Gleichungssystem löst.