• No results found

2.6. Vektorer och matriser i Fortran

N/A
N/A
Protected

Academic year: 2021

Share "2.6. Vektorer och matriser i Fortran"

Copied!
21
0
0

Loading.... (view fulltext now)

Full text

(1)

forts. p˚a f¨oreg˚aende f¨orel¨asning (EXIT och CYCLE):

I Fortran 90 ¨ar det inte n¨odv¨andigt att anv¨anda en DO-variabel, emedan det ocks˚a finns ett annat s¨att att konstruera en DO-slinga, d¨ar en ny sats, EXIT, anv¨ands f¨or att ¨overflytta kontrollen till den sats som f¨oljer n¨armast efter END DO-satsen. N¨ar EXIT-satsen utf¨ors, kommer s˚aledes alla programsatser som befinner sig mellan EXIT och END DO att hoppas ¨over. EXIT-satsen m˚aste d¨arf¨or f¨oreg˚as av en IF-sats, t.ex. p˚a f¨oljande s¨att:

DO ...

IF (x < eps) EXIT ...

END DO

! programmet forts¨atter h¨ar, efter EXIT ...

Detta s¨att att konstruera en DO-slinga ¨ar dock riskabelt, eftersom programmet kan r˚aka ut f¨or en o¨andlig slinga, som man inte kan stoppa, om inte programmet (eller datorn) stoppas med v˚ald! F¨or s¨akerhets skull

¨ar det d¨arf¨or alltid b¨ast att anv¨anda EXIT-satsen i kombination med en DO-variabel, vars gr¨anser ¨ar angivna, t.ex. p˚a f¨oljande s¨att:

(2)

DO ind = min, max ...

IF (x < eps) EXIT ...

END DO

! programmet forts¨atter h¨ar, efter EXIT, eller d˚a maximiantalet

! iterationer har n˚atts

En f¨ordel med denna metod ¨ar, att man dessutom k¨anner till hur m˚anga iterationer som blivit utf¨orda, eftersom ind kommer att f˚a v¨ardet max+1 om slingan g˚atts igenom i sin helhet.

Som ett exempel p˚a anv¨andningen av EXIT ger vi f¨oljande program, som r¨aknar ut f¨ordelningen av positiva och negativa tal bland slumpm¨assigt valda heltal. R¨akningen avbryts antingen d˚a talet 0 n˚as, eller d˚a ett p˚a f¨orhand givet antal slumptal ber¨aknats.

PROGRAM rndtest IMPLICIT NONE

CHARACTER (LEN=7) :: kod REAL :: rnd

INTEGER :: num, nr_pos, nr_neg, nr_max, i nr_pos = 0

nr_neg = 0

PRINT *, "Ge antal slumptal:"

READ *, nr_max

(3)

DO i=1, nr_max

CALL RANDOM_NUMBER(rnd) num = nr_max*(rnd-0.5) kod = " "

IF (num < 0) THEN kod = "negativ"

ELSE IF (num > 0) THEN kod = "positiv"

END IF

SELECT CASE (kod) CASE ("positiv")

nr_pos = nr_pos + 1 CASE ("negativ")

nr_neg= nr_neg + 1 CASE DEFAULT

EXIT END SELECT END DO

PRINT *, "Antal positiva tal:", nr_pos PRINT *, "Antal negativa tal:", nr_neg STOP

END PROGRAM rndtest

(4)

Som ett annat exempel skall vi studera ett program, som approximerar den best¨amda integralen av kurvan y = f (x) ¨over intervallet [a, b] med trapetsregeln. Som vi tidigare sett, s˚a kan vi indela intervallet [a, b]

i n lika stora delintervall, och konstruera n trapets, vilkas baser sammanfaller med delintervallen, och de andra h¨ornpunkterna ligger p˚a kurvan. Trapetsens sammanlagda yta kan d˚a uttryckas

In = h f (a)

2 + f (a + h) + f (a + 2h) + . . . + f (b − h) + f (b) 2

 ,

d¨ar h = (b − a)/n ¨ar l¨angden av varje delintervall. I programmet Integral ber¨aknas denna summa genom att man f¨orst bildar 12(f (a) + f (b)) och sedan adderar till de ¨ovriga ordinatorna. Programmet har till¨ampats p˚a funktionen y = sin(x). N¨ar skillnaden mellan approximationerna till integralen blivit tillr¨ackligt liten, avbryts r¨akningen med EXIT–kommandot.

PROGRAM Integral

! Ber¨akna en approximation till en integral

! genom att dela ett intervall [a,b] i N delintervall

! och till¨ampa trapetsregeln.

IMPLICIT NONE

REAL :: A,B,H,SUM,SUM1,EPS INTEGER :: N,I

EPS = 1.E-6

PRINT *, "Ange A, B:"

READ *, A,B

(5)

SUM1=0

DO N=100,10000,100 H = (B-A)/N

SUM = 0.5*(SIN(A)+SIN(B)) DO I = 1, N-1

SUM = SUM + SIN(A+I*H) END DO

SUM = H*SUM

PRINT *, ’N = ’,N, ’ SUM = ’,SUM IF (ABS(SUM-SUM1) < EPS) EXIT SUM1 = SUM

END DO STOP

END PROGRAM Integral

Som ett exempel skall vi ber¨akna integralen ¨over intervallet [0, π]:

Ange A, B: 0,3.14159

N = 100 SUM = 1.99984

N = 200 SUM = 1.99996

N = 300 SUM = 1.99998

N = 400 SUM = 1.99999

N = 500 SUM = 1.99999

N = 600 SUM = 2.00000

N = 700 SUM = 2.00000

(6)

I vissa fall har man behov av att flytta programkontrollen tillbaka till b¨orjan av DO-slingan, utan att behandla alla satserna i slingan. Detta ¨ar m¨ojligt med satsen CYCLE, vars anv¨andning l¨attast framg˚ar av ett exempel, d¨ar vi antar att programmet l¨aser in data i en DO-slinga, som behandlas vidare i en CASE-konstruktion, som ing˚ar i slingan. Vissa data antages leda till att programmet avslutas, medan felskrivningar leder till att data m˚aste l¨asas in p˚a nytt.

...

! Data har tagit slut CASE ("END")

EXIT

! Fel i inmatade data CASE DEFAULT

PRINT *, "Felaktiv kod, ge in data p˚a nytt"

CYCLE END SELECT ...

Den allm¨anna regeln ¨ar, att EXIT-satsen ¨overf¨or kontrollen till END DO-satsen i den innersta slingan d¨ar EXIT-satsen f¨orekommer. Detta g¨aller ocks˚a ifall vi har flere DO-slingor innanf¨or varandra, som i nedanst˚aende exempel:

(7)

DO ...

DO ...

DO ...

EXIT ...

END DO ...

END DO ...

END DO

Motsvarande regel g¨aller CYCLE-satsen. Ibland vill man kanske avsluta flera ¨an den innersta slingan. Detta

¨ar m¨ojligt i Fortran 90, emedan man kan ge namn ˚at DO-slingorna. Efter EXIT och CYCLE kan man sedan skriva namnet p˚a den slinga, ur vilken man hoppar ut, resp. upprepar. Detta visas tydligare av ett exempel:

(8)

slinga_1: DO ...

slinga_2: DO ...

SELECT CASE (ind) CASE (10)

EXIT slinga_2 CASE (20)

EXIT slinga_1 CASE (30)

CYCLE slinga_2 END SELECT

...

END DO slinga_2 ...

END DO slinga_1

I ovanst˚aende exempel ¨ar referenserna till slinga_2 egentligen ¨overfl¨odiga, men de g¨or programmet l¨attare att begripa. P˚a ett liknande s¨att kan man ocks˚a namnge IF- och CASE-konstruktioner f¨or att g¨ora kom- plicerade strukturer l¨attare att f¨orst˚a.

(9)

Som vi kan se, kan man genom att anv¨anda de nya formerna av DO-satserna ganska ofta undvika GO TO satser (som kr¨aver satsnummer och ¨ar det klassiska s¨attet att ¨andra programfl¨odets riktning), ¨aven om det fortfarande finns fall, d˚a programstrukturen d¨arigenom kan f¨orenklas.

En variant av DO-satsen, som redan ingick i m˚anga FORTRAN 77 kompilatorer som en utvidgning av spr˚aket, ¨ar DO WHILE (logiskt uttryck). Ett exempel ¨ar

fac = 1

DO WHILE (n>0) fac = fac*n n = n-1

END DO

DO WHILE–konstruktioner kan emellertid ers¨attas av DO

IF .NOT. (logiskt uttryck) EXIT ...

som har den f¨ordelen, att EXIT-satsen kan placeras var som helst i DO-slingan.

(10)

Vi skall ¨annu studera ett annat exempel, n¨amligen Euklides algoritm f¨or att ber¨akna st¨orsta gemensamma delaren, eller divisorn (sgd) f¨or tv˚a tal, som kan uttryckas s˚ah¨ar: 1) Dividera m med n och beteckna resten vid divisionen med r. 2) Om r = 0 s˚a ¨ar ber¨akningen klar och resultatet finns i m. 3) S¨att annars m till n och n till r. Upprepa fr˚an punkt 1.

PROGRAM euklides

! Programmet ber¨aknar sgd f¨or tv˚a heltal m och n med Euklides algoritm.

IMPLICIT NONE

INTEGER :: m, n, r

PRINT *, ’Ange tv˚a positiva heltal m och n:’

READ *, m,n

PRINT *, ’Talen ¨ar ’,m, ’ och ’,n IF (m>0 .AND. n>0) THEN

DO WHILE (n /= 0) r = MOD(m,n) m = n

n = r END DO

PRINT *, ’St¨orsta gemensamma divisorn: ’,m ELSE

PRINT *, ’Talen ¨ar inte positiva!’

END IF

END PROGRAM euklides

(11)

2.6. Vektorer och matriser i Fortran

N¨ar man behandlar data i fysiken, anv¨ander man ofta vektor- eller matrisframst¨allning. F¨or att behandla dessa ordnade strukturer, har Fortran indicerade variabler, som deklareras genom att man specificerar en maximistorlek, eller dimension i samband med variabelnamnet. Man kan ocks˚a deklarera vektorer och matriser med s.k. dimensionsattribut. H¨ar nedan visas n˚agra exempel p˚a s˚adana typdeklarationer:

REAL, DIMENSION(100) :: a, d(20) REAL :: b(100,100)

REAL, DIMENSION(10,10) :: mat INTEGER, DIMENSION(2,2) :: mat_2

Dessa deklarationer reserverar s˚alunda utrymme f¨or den reella vektorn a med 100 element, vektorn d med 20 element, en reell matris b med 100 rader och 100 kolumner, en reell matris mat med 10 rader och 10 kolumner samt en heltalsmatris mat_2 med tv˚a rader och tv˚a kolumner.

Dimensionstalet ¨ar ¨ovre gr¨ans f¨or index, och dess undre gr¨ans ¨ar 1 (utg˚angsantagande). Ibland kan det vara bekv¨amt att ange b˚ade en undre och en ¨ovre gr¨ans f¨or index skilt f¨or sig. Dessa gr¨anser skiljes i s˚adant fall ˚at av ett kolon: undre gr¨ans : ¨ovre gr¨ans. En vektor vars f¨orsta element ¨ar 0, och sista element 10, kan s˚alunda dimensioneras med instruktionen REAL, DIMENSION (0:10) :: a.

(12)

Indexgr¨anserna f˚ar ocks˚a vara aritmetiska uttryck, f¨orutsatt att v¨ardet av dessa uttryck ¨ar best¨amda. Efter- som dimensionssatserna och typdeklarationerna m˚aste st˚a f¨orst i programmet, kan dessa uttryck inte in- neh˚alla variabler. D¨aremot kan de inneh˚alla parameterkonstanter, som definieras med hj¨alp av PARAMETER- attributet, t. ex.

INTEGER, PARAMETER :: NMAX=1000 REAL, DIMENSION(NMAX/2) :: a

PARAMETER-attributet anger att ifr˚agavarande parameter har ett konstant v¨arde, som inte kan ¨andras under programexekveringen. Som av exemplet framg˚ar, m˚aste en parameterkonstant ocks˚a deklareras innan den anv¨ands f¨orsta g˚angen.

Man refererar till en matris normalt genom att ange ett element, t. ex. b(1,1), som betecknar matrisen b:s element i f¨orsta raden och f¨orsta kolumnen. Vektorn a och matrisen b kan vi generera p˚a f¨oljande s¨att:

DO i=1,100

a(i) = 110./(i+10.) END DO

DO i=1,10 DO j=1,10

b(i,j) = a((i-1)*10 + j) END DO

END DO

(13)

Som ett annat exempel kan vi v¨alja ett program som l¨aser in en talr¨acka och s¨oker upp det st¨orsta talet:

PROGRAM maxtal IMPLICIT NONE

REAL :: a(100),amax INTEGER :: i,n

! L¨as in talen:

PRINT *, ’Antalet tal (h¨ogst 100):’

READ *, n

PRINT *, ’Skriv talen:’

DO i=1,n

READ *, a(i) END DO

! Best¨am maximum:

amax = a(1) DO i=2,n

IF (a(i) > amax) amax = a(i) END DO

! Skriv ut resultatet:

PRINT *, ’Maximum ¨ar ’, amax STOP

END PROGRAM maxtal

(14)

Ofta har man behov att ordna elementen i en talr¨acka i storleksordning, dvs sortera dem. Detta kan man g¨ora s˚a, att man f¨orst s¨oker upp det st¨orsta talet i r¨ackan, och l˚ater det byta plats med det f¨orsta talet.

D¨arp˚a best¨ammer man maximum av talen i resten av talr¨ackan, och upprepar proceduren, ¨anda tills alla talen har ordnats i storleksordning. H¨arnedan visas ett enkelt sorteringsprogram enligt denna metod.

PROGRAM sort IMPLICIT NONE

REAL :: a(100), amax INTEGER :: i, j, n, imax

! Skriv in talen:

PRINT *, ’Ange antalet tal:’

READ *, n

PRINT *, ’Skriv talen:’

DO i=1,n

READ *, a(i) END DO

! Sortera talen:

DO i=1,n-1 amax = a(i) imax = i

! S¨ok upp maximum av talen i den osorterade

! delen av talr¨ackan:

(15)

DO j=i+1,n

IF (a(j) > amax) THEN amax = a(j)

imax = j END IF

END DO

a(imax) = a(i) a(i) = amax END DO

! Skriv ut resultatet:

DO i =1,n,5

PRINT *, (a(j), j= i,i+4) END DO

STOP

END PROGRAM sort

Resultatet skrivs ut med fem tal per rad. Observera, hur man kan kan skriva ut flera indicerade variabler p˚a samma rad ( s.k. implicerat DO).

(16)

Vi skall ¨annu visa hur man kan ber¨akna och skriva ut Pascals triangel. Talen, dvs binomialkoefficienterna, som triangeln uppbyggs av, kan man l¨att generera med hj¨alp av rekursionsformeln nk = n−1k−1 + n−1k .

PROGRAM Pascal IMPLICIT NONE

INTEGER :: c(0:9), k, n DO k=1,9

c(k)=0 END DO c(0) = 1 DO n=0,9

DO k=n,1,-1

c(k) = c(k-1) + c(k) END DO

PRINT *, (c(k), k=0,n) END DO

STOP

END PROGRAM Pascal

(17)

Observera, att elementen av en tv˚adimensionell matris alltid lagras kolumnvis, s˚a att t.ex. elementen f¨or en 3 × 2-matris B,

B =

B11 B12 B21 B22 B31 B32

 som deklarerats p˚a f¨oljande s¨att:

REAL, DIMENSION (1:3, 1:2) :: B

lagras i minnet i ordningsf¨oljden B(1,1), B(2,1), B(3,1), B(1,2), B(2,2), B(3,2) (dvs f¨orsta index varierar alltid snabbast).

F¨or att tilldela v¨arden ˚at elementen av en matris och en vektor, kan man anv¨anda tilldelningssatser, s˚asom ovan visats. Ibland vill man initialisera sina variabler, dvs ge utg˚angsv¨arden ˚at elementen. I Fortran 90 finns det ett beh¨andigt s¨att att konstruera matriser och vektorer, utg˚aende fr˚an en lista av matriselementen. Om man t.ex. vill konstruera en heltalsvektor med elementen [1 2 3 4 5], kan man anv¨anda satserna

INTEGER, DIMENSION (1:5) :: vec vec = (/1,2,3,4,5/)

Om man vill konstruera matriser av h¨ogre dimension, m˚aste man anv¨anda den inbyggda funktionen RESHAPE. Det enklaste s¨attet att anv¨anda den ¨ar genom anropet

(18)

Matris = RESHAPE(vec1, vec2)

d¨ar vec1 ¨ar en vektor, som inneh˚aller matriselementen, och vec2 ¨ar en vektor, som anger matrisens form.

Som ett exempel skall vi konstruera matrisen

A =

1 2 3

4 5 6

7 8 9

 vilket g˚ar till p˚a f¨oljande s¨att:

INTEGER :: vec(9)

INTEGER, DIMENSION(3,3) :: A vec = (/1,4,7,2,5,8,3,6,9/) A = RESHAPE(vec, (/3,3/))

Om matriselementen enkelt kan alstras genom n˚agon regel, s˚a kan man ibland ist¨allet f¨or att r¨akna upp alla elementen, generera dem via implicerat DO. I exemplet ovan kan matriselementen t.ex. genereras med instruktionen

(19)

vec = (/((i+j, j=1,7,3),i=0,2)/)

Den kolumnvisa lagringen leder ibland till missf¨orst˚and. L˚at oss t.ex. se p˚a den tv˚adimensionella matrisen A i ovanst˚aende exempel. Om man skriver ut den med PRINT-satsen PRINT *, A kommer kolumnerna att skrivas ut efter varandra, vilket inte var meningen. F¨or att f˚a matrisen utskriven radvis, m˚aste man anv¨anda en DO-slinga:

DO i = 1,3

PRINT *, (A(i,j), j=1,3) END DO

Ist¨allet f¨or implicerat DO: PRINT *, (A(i,j), j=1,3), skulle man i detta fall ocks˚a kunna skriva n˚agot kortare PRINT *, A(i,:), eftersom kolonet ers¨atter alla kolumner, som reserverats enligt matrisens typdeklaration. Observera, att A(:,:) ¨ar ekvivalent med A.

Man brukar ofta transponera matriser, dvs l˚ata kolumner och rader byta plats. Att transponera en matris i FORTRAN ¨ar inte helt trivialt, eftersom man b¨or vara f¨orsiktig, s˚a att man inte skriver ¨over matrisen.

(20)

En enkel metod att transponera en N × N matris visas av f¨oljande del av ett program:

DO I=1,N-1 DO J=I+1,N

TMP = A(I,J) A(I,J) = A(J,I) A(J,I) = TMP END DO

END DO

I detta program har ¨aven beaktats att man inte beh¨over transponera elementen p˚a matrisens diagonal. DO- variablernas gr¨anser har valts s˚a, att man inte i on¨odan upprepar samma operationer p˚a matrisens element.

F¨or att g¨ora det enklare att transponera matriser finns det i Fortran 90 en s¨arskild inbyggd funktion, som f¨orklaras b¨ast med ett exempel: B = TRANSPOSE(A).

Ett annat exempel p˚a matriser ¨ar ett program, som ber¨aknar produkten av en m × p-matris A och en p × n-matris B enligt formeln

Cij =

p

X

k=1

AikBkj. Matrisen C ¨ar s˚aledes en m × n-matris.

(21)

En slinga som utf¨or denna multiplikation kan skrivas p˚a f¨oljande s¨att: (Observera nollst¨allningen av variabeln sum, d¨ar produkterna lagras!)

DO i=1,m DO j=1,n

sum = 0.

DO k=1,p

sum = sum + A(i,k)*B(k,j) END DO

C(i,j) = sum END DO

END DO

F¨or att f¨orenkla matrismultiplikationen, finns det i Fortran 90 en s¨arskild funktion, n¨amligen MATMUL. Med hj¨alp av denna funktion kan matriserna i exemplet ovan multipliceras p˚a f¨oljande s¨att: C = MATMUL(A,B).

References

Related documents

Detta g¨aller alla tal vars dyadiska utveckling ¨ar ¨andlig; man beh¨over inte kasta fler kast ¨an vad som anges av den position d¨ar sista ettan finns i utvecklingen.. Det betyder

Till exempel fick jag inte med n˚ agot Ljus- och Optikland i f¨ orsta f¨ ors¨ oket, och pilen mot Kosmologi, som ligger utanf¨ or den h¨ ar kartan, borde peka mer upp˚ at,

Det ¨ ar en mots¨ agelse till att vi f˚ ar stryka alla gemensamma faktorer och d¨ arf¨ or ¨ ar x irrationellt.. (a) Skissa grafen av den trigonometriska

L˚ at y(t) vara andelen av populationen som ¨ar smittad efter tiden t dygn, r¨aknad fr˚ an uppt¨ack- ten... Observera att ¨amnets koncentration ¨ar samma som m¨angden av

¨ overf¨ oras till en annan dator, och assemblerprogrammen dessutom oftast ¨ar l˚ anga och invecklade, ins˚ ag man redan tidigt vikten av att uppfinna h¨ ogre programmeringsspr˚

F¨ or att formatera teckenstr¨angar anv¨ands i Fortran ett s¨arskilt format, som allm¨ant kan uttryckas kAn, d¨ar k anger antalet upprepningar, och n ¨ar antalet tecken i

Det ¨ar d¨arf¨ or b¨attre om man kan kontrollera anv¨andningen av komponenterna av en h¨arledd typ i programmet, s˚ a att de antingen ¨ar fritt tillg¨angliga i hela

Rutinen som anv¨ands f¨ or att definiera operatorn, kan ha antingen ett eller tv˚ a argument, men eftersom funktionen normalt definieras i samma modul som inneh˚