1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
|
c subroutine invers (a, max, n, iflag)
c
c Although it seems counter-intuitive, the tests that I have run
c so far suggest that the 180 x 180 matrices that NSTAR needs can
c be inverted with sufficient accuracy if the elements are REAL*4
c rather than REAL*8.
c
c Arguments
c
c a (input/output) is a square matrix of dimension N. The inverse
c of the input matrix A is returned in A.
c
c max (input) is the size assigned to the matrix A in the calling
c routine. It is needed for the dimension statement below.
c
c iflag (output) is an error flag. iflag = 1 if the matrix could not
c be inverted; iflag = 0 if it could.
c
subroutine invers (a, max, n, iflag)
c
implicit none
integer max, n, iflag
real a(max,max)
integer i, j, k
c
iflag = 0
i = 1
300 if (a(i,i) .eq. 0.0e0) go to 9100
a(i,i) = 1.0e0 / a(i,i)
j = 1
301 if (j .eq. i) go to 304
a(j,i) = -a(j,i) * a(i,i)
k = 1
302 if (k .eq. i) go to 303
a(j,k) = a(j,k) + a(j,i) * a(i,k)
303 if (k .eq. n) go to 304
k = k + 1
go to 302
304 if (j .eq. n) go to 305
j = j + 1
go to 301
305 k = 1
306 if (k .eq. i) go to 307
a(i,k) = a(i,k) * a(i,i)
307 if (k .eq. n) go to 308
k = k + 1
go to 306
308 if (i .eq. n) return
i = i+1
go to 300
c
c Error: zero on the diagonal.
c
9100 iflag = 1
return
c
end
c
c
c
c subroutine dinvers (a, max, n, iflag)
c
c Arguments
c
c a (input/output) is a square matrix of dimension N. The inverse
c of the input matrix A is returned in A.
c
c max (input) is the size assigned to the matrix A in the calling
c routine. It's needed for the dimension statement below.
c
c iflag (output) is an error flag. iflag = 1 if the matrix could not
c be inverted; iflag = 0 if it could.
c
subroutine dinvers (a, max, n, iflag)
c
implicit none
integer max, n, iflag
double precision a(max,max)
integer i, j, k
c
iflag = 0
i = 1
300 if (a(i,i) .eq. 0.0e0) go to 9100
a(i,i) = 1.0e0 / a(i,i)
j = 1
301 if (j .eq. i) go to 304
a(j,i) = -a(j,i) * a(i,i)
k = 1
302 if (k .EQ. i) go to 303
a(j,k) = a(j,k) + a(j,i) * a(i,k)
303 if (k .eq. n) go to 304
k = k + 1
go to 302
304 if (j .eq. n) go to 305
j = j + 1
go to 301
305 k = 1
306 if (k .eq. i) go to 307
a(i,k) = a(i,k) * a(i,i)
307 if (k .eq. n) go to 308
k = k + 1
go to 306
308 if (i .eq. n) return
i = i+1
go to 300
c
c Error: zero on the diagonal.
c
9100 iflag = 1
return
c
end
|
