mirror of
https://github.com/Cantera/cantera.git
synced 2025-02-25 18:55:29 -06:00
136 lines
3.7 KiB
Fortran
136 lines
3.7 KiB
Fortran
subroutine dgbsl(abd,lda,n,ml,mu,ipvt,b,job)
|
|
integer lda,n,ml,mu,ipvt(1),job
|
|
double precision abd(lda,1),b(1)
|
|
c
|
|
c dgbsl solves the double precision band system
|
|
c a * x = b or trans(a) * x = b
|
|
c using the factors computed by dgbco or dgbfa.
|
|
c
|
|
c on entry
|
|
c
|
|
c abd double precision(lda, n)
|
|
c the output from dgbco or dgbfa.
|
|
c
|
|
c lda integer
|
|
c the leading dimension of the array abd .
|
|
c
|
|
c n integer
|
|
c the order of the original matrix.
|
|
c
|
|
c ml integer
|
|
c number of diagonals below the main diagonal.
|
|
c
|
|
c mu integer
|
|
c number of diagonals above the main diagonal.
|
|
c
|
|
c ipvt integer(n)
|
|
c the pivot vector from dgbco or dgbfa.
|
|
c
|
|
c b double precision(n)
|
|
c the right hand side vector.
|
|
c
|
|
c job integer
|
|
c = 0 to solve a*x = b ,
|
|
c = nonzero to solve trans(a)*x = b , where
|
|
c trans(a) is the transpose.
|
|
c
|
|
c on return
|
|
c
|
|
c b the solution vector x .
|
|
c
|
|
c error condition
|
|
c
|
|
c a division by zero will occur if the input factor contains a
|
|
c zero on the diagonal. technically this indicates singularity
|
|
c but it is often caused by improper arguments or improper
|
|
c setting of lda . it will not occur if the subroutines are
|
|
c called correctly and if dgbco has set rcond .gt. 0.0
|
|
c or dgbfa has set info .eq. 0 .
|
|
c
|
|
c to compute inverse(a) * c where c is a matrix
|
|
c with p columns
|
|
c call dgbco(abd,lda,n,ml,mu,ipvt,rcond,z)
|
|
c if (rcond is too small) go to ...
|
|
c do 10 j = 1, p
|
|
c call dgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0)
|
|
c 10 continue
|
|
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,ddot
|
|
c fortran min0
|
|
c
|
|
c internal variables
|
|
c
|
|
double precision ddot,t
|
|
integer k,kb,l,la,lb,lm,m,nm1
|
|
c
|
|
m = mu + ml + 1
|
|
nm1 = n - 1
|
|
if (job .ne. 0) go to 50
|
|
c
|
|
c job = 0 , solve a * x = b
|
|
c first solve l*y = b
|
|
c
|
|
if (ml .eq. 0) go to 30
|
|
if (nm1 .lt. 1) go to 30
|
|
do 20 k = 1, nm1
|
|
lm = min0(ml,n-k)
|
|
l = ipvt(k)
|
|
t = b(l)
|
|
if (l .eq. k) go to 10
|
|
b(l) = b(k)
|
|
b(k) = t
|
|
10 continue
|
|
call daxpy(lm,t,abd(m+1,k),1,b(k+1),1)
|
|
20 continue
|
|
30 continue
|
|
c
|
|
c now solve u*x = y
|
|
c
|
|
do 40 kb = 1, n
|
|
k = n + 1 - kb
|
|
b(k) = b(k)/abd(m,k)
|
|
lm = min0(k,m) - 1
|
|
la = m - lm
|
|
lb = k - lm
|
|
t = -b(k)
|
|
call daxpy(lm,t,abd(la,k),1,b(lb),1)
|
|
40 continue
|
|
go to 100
|
|
50 continue
|
|
c
|
|
c job = nonzero, solve trans(a) * x = b
|
|
c first solve trans(u)*y = b
|
|
c
|
|
do 60 k = 1, n
|
|
lm = min0(k,m) - 1
|
|
la = m - lm
|
|
lb = k - lm
|
|
t = ddot(lm,abd(la,k),1,b(lb),1)
|
|
b(k) = (b(k) - t)/abd(m,k)
|
|
60 continue
|
|
c
|
|
c now solve trans(l)*x = y
|
|
c
|
|
if (ml .eq. 0) go to 90
|
|
if (nm1 .lt. 1) go to 90
|
|
do 80 kb = 1, nm1
|
|
k = n - kb
|
|
lm = min0(ml,n-k)
|
|
b(k) = b(k) + ddot(lm,abd(m+1,k),1,b(k+1),1)
|
|
l = ipvt(k)
|
|
if (l .eq. k) go to 70
|
|
t = b(l)
|
|
b(l) = b(k)
|
|
b(k) = t
|
|
70 continue
|
|
80 continue
|
|
90 continue
|
|
100 continue
|
|
return
|
|
end
|