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

Annotation of /trunk/linpack/dgbfa.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (hide annotations) (download)
Fri Oct 29 20:54:12 2004 UTC (15 years, 9 months ago) by aw0a
File size: 5418 byte(s)
Setting up web subdirectory in repository
1 aw0a 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