/[ascend]/trunk/linpack/dgbfa.f
ViewVC logotype

Contents of /trunk/linpack/dgbfa.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download)
Fri Oct 29 20:54:12 2004 UTC (20 years ago) by aw0a
File size: 5418 byte(s)
Setting up web subdirectory in repository
1 C dgbfa.f
2 C is freely available from netlib. It is not subject to any GNU License
3 C set by the authors of the ASCEND math programming system.
4 C $Date: 1996/04/30 18:17:11 $ $Revision: 1.1.1.1 $
5 C
6 subroutine dgbfa(abd,lda,n,ml,mu,ipvt,info)
7 integer lda,n,ml,mu,ipvt(1),info
8 double precision abd(lda,1)
9 c
10 c dgbfa factors a double precision band matrix by elimination.
11 c
12 c dgbfa is usually called by dgbco, but it can be called
13 c directly with a saving in time if rcond is not needed.
14 c
15 c on entry
16 c
17 c abd double precision(lda, n)
18 c contains the matrix in band storage. the columns
19 c of the matrix are stored in the columns of abd and
20 c the diagonals of the matrix are stored in rows
21 c ml+1 through 2*ml+mu+1 of abd .
22 c see the comments below for details.
23 c
24 c lda integer
25 c the leading dimension of the array abd .
26 c lda must be .ge. 2*ml + mu + 1 .
27 c
28 c n integer
29 c the order of the original matrix.
30 c
31 c ml integer
32 c number of diagonals below the main diagonal.
33 c 0 .le. ml .lt. n .
34 c
35 c mu integer
36 c number of diagonals above the main diagonal.
37 c 0 .le. mu .lt. n .
38 c more efficient if ml .le. mu .
39 c on return
40 c
41 c abd an upper triangular matrix in band storage and
42 c the multipliers which were used to obtain it.
43 c the factorization can be written a = l*u where
44 c l is a product of permutation and unit lower
45 c triangular matrices and u is upper triangular.
46 c
47 c ipvt integer(n)
48 c an integer vector of pivot indices.
49 c
50 c info integer
51 c = 0 normal value.
52 c = k if u(k,k) .eq. 0.0 . this is not an error
53 c condition for this subroutine, but it does
54 c indicate that dgbsl will divide by zero if
55 c called. use rcond in dgbco for a reliable
56 c indication of singularity.
57 c
58 c band storage
59 c
60 c if a is a band matrix, the following program segment
61 c will set up the input.
62 c
63 c ml = (band width below the diagonal)
64 c mu = (band width above the diagonal)
65 c m = ml + mu + 1
66 c do 20 j = 1, n
67 c i1 = max0(1, j-mu)
68 c i2 = min0(n, j+ml)
69 c do 10 i = i1, i2
70 c k = i - j + m
71 c abd(k,j) = a(i,j)
72 c 10 continue
73 c 20 continue
74 c
75 c this uses rows ml+1 through 2*ml+mu+1 of abd .
76 c in addition, the first ml rows in abd are used for
77 c elements generated during the triangularization.
78 c the total number of rows needed in abd is 2*ml+mu+1 .
79 c the ml+mu by ml+mu upper left triangle and the
80 c ml by ml lower right triangle are not referenced.
81 c
82 c linpack. this version dated 08/14/78 .
83 c cleve moler, university of new mexico, argonne national lab.
84 c
85 c subroutines and functions
86 c
87 c blas daxpy,dscal,idamax
88 c fortran max0,min0
89 c
90 c internal variables
91 c
92 double precision t
93 integer i,idamax,i0,j,ju,jz,j0,j1,k,kp1,l,lm,m,mm,nm1
94 c
95 c
96 m = ml + mu + 1
97 info = 0
98 c
99 c zero initial fill-in columns
100 c
101 j0 = mu + 2
102 j1 = min0(n,m) - 1
103 if (j1 .lt. j0) go to 30
104 do 20 jz = j0, j1
105 i0 = m + 1 - jz
106 do 10 i = i0, ml
107 abd(i,jz) = 0.0d0
108 10 continue
109 20 continue
110 30 continue
111 jz = j1
112 ju = 0
113 c
114 c gaussian elimination with partial pivoting
115 c
116 nm1 = n - 1
117 if (nm1 .lt. 1) go to 130
118 do 120 k = 1, nm1
119 kp1 = k + 1
120 c
121 c zero next fill-in column
122 c
123 jz = jz + 1
124 if (jz .gt. n) go to 50
125 if (ml .lt. 1) go to 50
126 do 40 i = 1, ml
127 abd(i,jz) = 0.0d0
128 40 continue
129 50 continue
130 c
131 c find l = pivot index
132 c
133 lm = min0(ml,n-k)
134 l = idamax(lm+1,abd(m,k),1) + m - 1
135 ipvt(k) = l + k - m
136 c
137 c zero pivot implies this column already triangularized
138 c
139 if (abd(l,k) .eq. 0.0d0) go to 100
140 c
141 c interchange if necessary
142 c
143 if (l .eq. m) go to 60
144 t = abd(l,k)
145 abd(l,k) = abd(m,k)
146 abd(m,k) = t
147 60 continue
148 c
149 c compute multipliers
150 c
151 t = -1.0d0/abd(m,k)
152 call dscal(lm,t,abd(m+1,k),1)
153 c
154 c row elimination with column indexing
155 c
156 ju = min0(max0(ju,mu+ipvt(k)),n)
157 mm = m
158 if (ju .lt. kp1) go to 90
159 do 80 j = kp1, ju
160 l = l - 1
161 mm = mm - 1
162 t = abd(l,j)
163 if (l .eq. mm) go to 70
164 abd(l,j) = abd(mm,j)
165 abd(mm,j) = t
166 70 continue
167 call daxpy(lm,t,abd(m+1,k),1,abd(mm+1,j),1)
168 80 continue
169 90 continue
170 go to 110
171 100 continue
172 info = k
173 110 continue
174 120 continue
175 130 continue
176 ipvt(n) = n
177 if (abd(m,n) .eq. 0.0d0) info = n
178 return
179 end

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22