Bemerkung: Das Array aa, das in Zeile 10 deklariert wird, ist ein sogenanntes automatisches Array, dessen Dimensionen tatsächlich erst zur Laufzeit durch den Wert des Parameters n bestimmt sind. Dadurch erübrigt sich das Übergeben eines Work-Arrays als Parameter oder - noch schlimmer - das Einführen eines künstlichen Parameters NMAX in GAUSS. Dies ist eine wichtige Neuerung in Fortran 90 und wird auch vom Convex-Fortran-Compiler unterstützt, ebenso wie viele andere Array-Erweiterungen von Fortran 90.
1 subroutine gauss (n, a, x, b )
2
3 implicit none
4 intrinsic abs
5
6 integer n
7 real*8 a(n,n), x(n), b(n)
8
9 integer i, j, k, ipivot
10 real*8 aa(n,n+1), fac, sto
11
12 real*8 eps
13 parameter(eps = 1.0d-10)
14
15 ! Umspeichern der Matrix a und der rechten Seite b
16 do i = 1, n
17 do j = 1, n
18 aa(i,j) = a(i,j)
19 enddo
20 aa(i,n+1) = b(i)
21 enddo
22
23 do i = 1, n-1
24
25 ! Pivotisierung
26 fac = 0.d0
27 ipivot = i
28 do j = i, n
29 if ( abs(aa(j,i)) .gt. fac ) then
30 fac = abs(aa(j,i))
31 ipivot = j
32 endif
33 enddo
34
35 if ( fac .lt. eps ) stop 'in gauss: Matrix singulaer'
36
37 if ( ipivot .ne. i ) then
38 do j = i, n+1
39 sto = aa(i,j)
40 aa(i,j) = aa(ipivot,j)
41 aa(ipivot,j) = sto
42 enddo
43 endif
44
45 ! Transformation der Matrix auf obere Dreiecksgestalt
46 do j = i+1, n
47 fac = aa(j,i) / aa(i,i)
48 do k = i+1, n+1
49 aa(j,k) = aa(j,k) - fac * aa(i,k)
50 enddo
51 enddo
52
53 enddo
54
55 ! Rueckwaertseinsetzen
56 x(n) = aa(n,n+1) / aa(n,n)
57 do i = n-1, 1, -1
58 x(i) = aa(i,n+1)
59 do j = n, i+1, -1
60 x(i) = x(i) - x(j) * aa(i,j)
61 enddo
62 x(i) = x(i) / aa(i,i)
63 enddo
64
65 return
66 end