mirror of
https://github.com/Cantera/cantera.git
synced 2025-02-25 18:55:29 -06:00
175 lines
5.1 KiB
Fortran
175 lines
5.1 KiB
Fortran
subroutine dgbfa(abd,lda,n,ml,mu,ipvt,info)
|
|
integer lda,n,ml,mu,ipvt(1),info
|
|
double precision abd(lda,1)
|
|
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 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 = max0(1, j-mu)
|
|
c i2 = min0(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 linpack. this version dated 08/14/78 .
|
|
c cleve moler, university of new mexico, argonne national lab.
|
|
c
|
|
c subroutines and functions
|
|
c
|
|
c blas daxpy,dscal,idamax
|
|
c fortran max0,min0
|
|
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 = min0(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.0d0
|
|
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.0d0
|
|
40 continue
|
|
50 continue
|
|
c
|
|
c find l = pivot index
|
|
c
|
|
lm = min0(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.0d0) 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.0d0/abd(m,k)
|
|
call dscal(lm,t,abd(m+1,k),1)
|
|
c
|
|
c row elimination with column indexing
|
|
c
|
|
ju = min0(max0(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.0d0) info = n
|
|
return
|
|
end
|