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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
|
subroutine dgbfa(abd,lda,n,ml,mu,ipvt,info)
integer lda,n,ml,mu,ipvt(*),info
double precision abd(lda,*)
c!purpose
c
c dgbfa factors a double precision band matrix by elimination.
c
c dgbfa is usually called by dgbco, but it can be called
c directly with a saving in time if rcond is not needed.
c
c!calling sequence
c
c subroutine dgbfa(abd,lda,n,ml,mu,ipvt,info)
c on entry
c
c abd double precision(lda, n)
c contains the matrix in band storage. the columns
c of the matrix are stored in the columns of abd and
c the diagonals of the matrix are stored in rows
c ml+1 through 2*ml+mu+1 of abd .
c see the comments below for details.
c
c lda integer
c the leading dimension of the array abd .
c lda must be .ge. 2*ml + mu + 1 .
c
c n integer
c the order of the original matrix.
c
c ml integer
c number of diagonals below the main diagonal.
c 0 .le. ml .lt. n .
c
c mu integer
c number of diagonals above the main diagonal.
c 0 .le. mu .lt. n .
c more efficient if ml .le. mu .
c on return
c
c abd an upper triangular matrix in band storage and
c the multipliers which were used to obtain it.
c the factorization can be written a = l*u where
c l is a product of permutation and unit lower
c triangular matrices and u is upper triangular.
c
c ipvt integer(n)
c an integer vector of pivot indices.
c
c info integer
c = 0 normal value.
c = k if u(k,k) .eq. 0.0 . this is not an error
c condition for this subroutine, but it does
c indicate that dgbsl will divide by zero if
c called. use rcond in dgbco for a reliable
c indication of singularity.
c
c band storage
c
c if a is a band matrix, the following program segment
c will set up the input.
c
c ml = (band width below the diagonal)
c mu = (band width above the diagonal)
c m = ml + mu + 1
c do 20 j = 1, n
c i1 = max(1, j-mu)
c i2 = min(n, j+ml)
c do 10 i = i1, i2
c k = i - j + m
c abd(k,j) = a(i,j)
c 10 continue
c 20 continue
c
c this uses rows ml+1 through 2*ml+mu+1 of abd .
c in addition, the first ml rows in abd are used for
c elements generated during the triangularization.
c the total number of rows needed in abd is 2*ml+mu+1 .
c the ml+mu by ml+mu upper left triangle and the
c ml by ml lower right triangle are not referenced.
c
c!originator
c linpack. this version dated 08/14/78 .
c cleve moler, university of new mexico, argonne national lab.
c
c!auxiliary routines
c
c blas daxpy,dscal,idamax
c fortran max,min
c
c!
c internal variables
c
double precision t
integer i,idamax,i0,j,ju,jz,j0,j1,k,kp1,l,lm,m,mm,nm1
c
c
m = ml + mu + 1
info = 0
c
c zero initial fill-in columns
c
j0 = mu + 2
j1 = min(n,m) - 1
if (j1 .lt. j0) go to 30
do 20 jz = j0, j1
i0 = m + 1 - jz
do 10 i = i0, ml
abd(i,jz) = 0.0d+0
10 continue
20 continue
30 continue
jz = j1
ju = 0
c
c gaussian elimination with partial pivoting
c
nm1 = n - 1
if (nm1 .lt. 1) go to 130
do 120 k = 1, nm1
kp1 = k + 1
c
c zero next fill-in column
c
jz = jz + 1
if (jz .gt. n) go to 50
if (ml .lt. 1) go to 50
do 40 i = 1, ml
abd(i,jz) = 0.0d+0
40 continue
50 continue
c
c find l = pivot index
c
lm = min(ml,n-k)
l = idamax(lm+1,abd(m,k),1) + m - 1
ipvt(k) = l + k - m
c
c zero pivot implies this column already triangularized
c
if (abd(l,k) .eq. 0.0d+0) go to 100
c
c interchange if necessary
c
if (l .eq. m) go to 60
t = abd(l,k)
abd(l,k) = abd(m,k)
abd(m,k) = t
60 continue
c
c compute multipliers
c
t = -1.0d+0/abd(m,k)
call dscal(lm,t,abd(m+1,k),1)
c
c row elimination with column indexing
c
ju = min(max(ju,mu+ipvt(k)),n)
mm = m
if (ju .lt. kp1) go to 90
do 80 j = kp1, ju
l = l - 1
mm = mm - 1
t = abd(l,j)
if (l .eq. mm) go to 70
abd(l,j) = abd(mm,j)
abd(mm,j) = t
70 continue
call daxpy(lm,t,abd(m+1,k),1,abd(mm+1,j),1)
80 continue
90 continue
go to 110
100 continue
info = k
110 continue
120 continue
130 continue
ipvt(n) = n
if (abd(m,n) .eq. 0.0d+0) info = n
return
end
|
