712 lines
21 KiB
Fortran
712 lines
21 KiB
Fortran
*DECK PCOEF
|
|
SUBROUTINE PCOEF (L, C, TC, A)
|
|
C***BEGIN PROLOGUE PCOEF
|
|
C***PURPOSE Convert the POLFIT coefficients to Taylor series form.
|
|
C***LIBRARY SLATEC
|
|
C***CATEGORY K1A1A2
|
|
C***TYPE SINGLE PRECISION (PCOEF-S, DPCOEF-D)
|
|
C***KEYWORDS CURVE FITTING, DATA FITTING, LEAST SQUARES, POLYNOMIAL FIT
|
|
C***AUTHOR Shampine, L. F., (SNLA)
|
|
C Davenport, S. M., (SNLA)
|
|
C***DESCRIPTION
|
|
C
|
|
C Written BY L. F. Shampine and S. M. Davenport.
|
|
C
|
|
C Abstract
|
|
C
|
|
C POLFIT computes the least squares polynomial fit of degree L as
|
|
C a sum of orthogonal polynomials. PCOEF changes this fit to its
|
|
C Taylor expansion about any point C , i.e. writes the polynomial
|
|
C as a sum of powers of (X-C). Taking C=0. gives the polynomial
|
|
C in powers of X, but a suitable non-zero C often leads to
|
|
C polynomials which are better scaled and more accurately evaluated.
|
|
C
|
|
C The parameters for PCOEF are
|
|
C
|
|
C INPUT --
|
|
C L - Indicates the degree of polynomial to be changed to
|
|
C its Taylor expansion. To obtain the Taylor
|
|
C coefficients in reverse order, input L as the
|
|
C negative of the degree desired. The absolute value
|
|
C of L must be less than or equal to NDEG, the highest
|
|
C degree polynomial fitted by POLFIT .
|
|
C C - The point about which the Taylor expansion is to be
|
|
C made.
|
|
C A - Work and output array containing values from last
|
|
C call to POLFIT .
|
|
C
|
|
C OUTPUT --
|
|
C TC - Vector containing the first LL+1 Taylor coefficients
|
|
C where LL=ABS(L). If L.GT.0 , the coefficients are
|
|
C in the usual Taylor series order, i.e.
|
|
C P(X) = TC(1) + TC(2)*(X-C) + ... + TC(N+1)*(X-C)**N
|
|
C If L .LT. 0, the coefficients are in reverse order,
|
|
C i.e.
|
|
C P(X) = TC(1)*(X-C)**N + ... + TC(N)*(X-C) + TC(N+1)
|
|
C
|
|
C***REFERENCES L. F. Shampine, S. M. Davenport and R. E. Huddleston,
|
|
C Curve fitting by polynomials in one variable, Report
|
|
C SLA-74-0270, Sandia Laboratories, June 1974.
|
|
C***ROUTINES CALLED PVALUE
|
|
C***REVISION HISTORY (YYMMDD)
|
|
C 740601 DATE WRITTEN
|
|
C 890531 Changed all specific intrinsics to generic. (WRB)
|
|
C 890531 REVISION DATE from Version 3.2
|
|
C 891214 Prologue converted to Version 4.0 format. (BAB)
|
|
C 920501 Reformatted the REFERENCES section. (WRB)
|
|
C***END PROLOGUE PCOEF
|
|
C
|
|
DIMENSION A(*), TC(*)
|
|
C***FIRST EXECUTABLE STATEMENT PCOEF
|
|
LL = ABS(L)
|
|
LLP1 = LL + 1
|
|
CALL PVALUE (LL,LL,C,TC(1),TC(2),A)
|
|
IF (LL .LT. 2) GO TO 2
|
|
FAC = 1.0
|
|
DO 1 I = 3,LLP1
|
|
FAC = FAC*(I-1)
|
|
1 TC(I) = TC(I)/FAC
|
|
2 IF (L .GE. 0) GO TO 4
|
|
NR = LLP1/2
|
|
LLP2 = LL + 2
|
|
DO 3 I = 1,NR
|
|
SAVE = TC(I)
|
|
NEW = LLP2 - I
|
|
TC(I) = TC(NEW)
|
|
3 TC(NEW) = SAVE
|
|
4 RETURN
|
|
END
|
|
c$$$
|
|
c$$$ subroutine dscal(n,da,dx,incx)
|
|
c$$$c
|
|
c$$$c scales a vector by a constant.
|
|
c$$$c uses unrolled loops for increment equal to one.
|
|
c$$$c jack dongarra, linpack, 3/11/78.
|
|
c$$$c modified 3/93 to return if incx .le. 0.
|
|
c$$$c
|
|
c$$$ double precision da,dx(1)
|
|
c$$$ integer i,incx,m,mp1,n,nincx
|
|
c$$$c
|
|
c$$$ if( n.le.0 .or. incx.le.0 )return
|
|
c$$$ if(incx.eq.1)go to 20
|
|
c$$$c
|
|
c$$$c code for increment not equal to 1
|
|
c$$$c
|
|
c$$$ nincx = n*incx
|
|
c$$$ do 10 i = 1,nincx,incx
|
|
c$$$ dx(i) = da*dx(i)
|
|
c$$$ 10 continue
|
|
c$$$ return
|
|
c$$$c
|
|
c$$$c code for increment equal to 1
|
|
c$$$c
|
|
c$$$c
|
|
c$$$c clean-up loop
|
|
c$$$c
|
|
c$$$ 20 m = mod(n,5)
|
|
c$$$ if( m .eq. 0 ) go to 40
|
|
c$$$ do 30 i = 1,m
|
|
c$$$ dx(i) = da*dx(i)
|
|
c$$$ 30 continue
|
|
c$$$ if( n .lt. 5 ) return
|
|
c$$$ 40 mp1 = m + 1
|
|
c$$$ do 50 i = mp1,n,5
|
|
c$$$ dx(i) = da*dx(i)
|
|
c$$$ dx(i + 1) = da*dx(i + 1)
|
|
c$$$ dx(i + 2) = da*dx(i + 2)
|
|
c$$$ dx(i + 3) = da*dx(i + 3)
|
|
c$$$ dx(i + 4) = da*dx(i + 4)
|
|
c$$$ 50 continue
|
|
c$$$ return
|
|
c$$$ end
|
|
|
|
subroutine dgbco(abd,lda,n,ml,mu,ipvt,rcond,z)
|
|
integer lda,n,ml,mu,ipvt(1)
|
|
double precision abd(lda,1),z(1)
|
|
double precision rcond
|
|
c
|
|
c dgbco factors a double precision band matrix by gaussian
|
|
c elimination and estimates the condition of the matrix.
|
|
c
|
|
c if rcond is not needed, dgbfa is slightly faster.
|
|
c to solve a*x = b , follow dgbco by dgbsl.
|
|
c to compute inverse(a)*c , follow dgbco by dgbsl.
|
|
c to compute determinant(a) , follow dgbco by dgbdi.
|
|
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
|
|
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 rcond double precision
|
|
c an estimate of the reciprocal condition of a .
|
|
c for the system a*x = b , relative perturbations
|
|
c in a and b of size epsilon may cause
|
|
c relative perturbations in x of size epsilon/rcond .
|
|
c if rcond is so small that the logical expression
|
|
c 1.0 + rcond .eq. 1.0
|
|
c is true, then a may be singular to working
|
|
c precision. in particular, rcond is zero if
|
|
c exact singularity is detected or the estimate
|
|
c underflows.
|
|
c
|
|
c z double precision(n)
|
|
c a work vector whose contents are usually unimportant.
|
|
c if a is close to a singular matrix, then z is
|
|
c an approximate null vector in the sense that
|
|
c norm(a*z) = rcond*norm(a)*norm(z) .
|
|
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 example.. if the original matrix is
|
|
c
|
|
c 11 12 13 0 0 0
|
|
c 21 22 23 24 0 0
|
|
c 0 32 33 34 35 0
|
|
c 0 0 43 44 45 46
|
|
c 0 0 0 54 55 56
|
|
c 0 0 0 0 65 66
|
|
c
|
|
c then n = 6, ml = 1, mu = 2, lda .ge. 5 and abd should contain
|
|
c
|
|
c * * * + + + , * = not used
|
|
c * * 13 24 35 46 , + = used for pivoting
|
|
c * 12 23 34 45 56
|
|
c 11 22 33 44 55 66
|
|
c 21 32 43 54 65 *
|
|
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 linpack dgbfa
|
|
c blas daxpy,ddot,dscal,dasum
|
|
c fortran dabs,dmax1,max0,min0,dsign
|
|
c
|
|
c internal variables
|
|
c
|
|
double precision ddot,ek,t,wk,wkm
|
|
double precision anorm,s,dasum,sm,ynorm
|
|
integer is,info,j,ju,k,kb,kp1,l,la,lm,lz,m,mm
|
|
c
|
|
c
|
|
c compute 1-norm of a
|
|
c
|
|
anorm = 0.0d0
|
|
l = ml + 1
|
|
is = l + mu
|
|
do 10 j = 1, n
|
|
anorm = dmax1(anorm,dasum(l,abd(is,j),1))
|
|
if (is .gt. ml + 1) is = is - 1
|
|
if (j .le. mu) l = l + 1
|
|
if (j .ge. n - ml) l = l - 1
|
|
10 continue
|
|
c
|
|
c factor
|
|
c
|
|
call dgbfa(abd,lda,n,ml,mu,ipvt,info)
|
|
c
|
|
c rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) .
|
|
c estimate = norm(z)/norm(y) where a*z = y and trans(a)*y = e .
|
|
c trans(a) is the transpose of a . the components of e are
|
|
c chosen to cause maximum local growth in the elements of w where
|
|
c trans(u)*w = e . the vectors are frequently rescaled to avoid
|
|
c overflow.
|
|
c
|
|
c solve trans(u)*w = e
|
|
c
|
|
ek = 1.0d0
|
|
do 20 j = 1, n
|
|
z(j) = 0.0d0
|
|
20 continue
|
|
m = ml + mu + 1
|
|
ju = 0
|
|
do 100 k = 1, n
|
|
if (z(k) .ne. 0.0d0) ek = dsign(ek,-z(k))
|
|
if (dabs(ek-z(k)) .le. dabs(abd(m,k))) go to 30
|
|
s = dabs(abd(m,k))/dabs(ek-z(k))
|
|
call dscal(n,s,z,1)
|
|
ek = s*ek
|
|
30 continue
|
|
wk = ek - z(k)
|
|
wkm = -ek - z(k)
|
|
s = dabs(wk)
|
|
sm = dabs(wkm)
|
|
if (abd(m,k) .eq. 0.0d0) go to 40
|
|
wk = wk/abd(m,k)
|
|
wkm = wkm/abd(m,k)
|
|
go to 50
|
|
40 continue
|
|
wk = 1.0d0
|
|
wkm = 1.0d0
|
|
50 continue
|
|
kp1 = k + 1
|
|
ju = min0(max0(ju,mu+ipvt(k)),n)
|
|
mm = m
|
|
if (kp1 .gt. ju) go to 90
|
|
do 60 j = kp1, ju
|
|
mm = mm - 1
|
|
sm = sm + dabs(z(j)+wkm*abd(mm,j))
|
|
z(j) = z(j) + wk*abd(mm,j)
|
|
s = s + dabs(z(j))
|
|
60 continue
|
|
if (s .ge. sm) go to 80
|
|
t = wkm - wk
|
|
wk = wkm
|
|
mm = m
|
|
do 70 j = kp1, ju
|
|
mm = mm - 1
|
|
z(j) = z(j) + t*abd(mm,j)
|
|
70 continue
|
|
80 continue
|
|
90 continue
|
|
z(k) = wk
|
|
100 continue
|
|
s = 1.0d0/dasum(n,z,1)
|
|
call dscal(n,s,z,1)
|
|
c
|
|
c solve trans(l)*y = w
|
|
c
|
|
do 120 kb = 1, n
|
|
k = n + 1 - kb
|
|
lm = min0(ml,n-k)
|
|
if (k .lt. n) z(k) = z(k) + ddot(lm,abd(m+1,k),1,z(k+1),1)
|
|
if (dabs(z(k)) .le. 1.0d0) go to 110
|
|
s = 1.0d0/dabs(z(k))
|
|
call dscal(n,s,z,1)
|
|
110 continue
|
|
l = ipvt(k)
|
|
t = z(l)
|
|
z(l) = z(k)
|
|
z(k) = t
|
|
120 continue
|
|
s = 1.0d0/dasum(n,z,1)
|
|
call dscal(n,s,z,1)
|
|
c
|
|
ynorm = 1.0d0
|
|
c
|
|
c solve l*v = y
|
|
c
|
|
do 140 k = 1, n
|
|
l = ipvt(k)
|
|
t = z(l)
|
|
z(l) = z(k)
|
|
z(k) = t
|
|
lm = min0(ml,n-k)
|
|
if (k .lt. n) call daxpy(lm,t,abd(m+1,k),1,z(k+1),1)
|
|
if (dabs(z(k)) .le. 1.0d0) go to 130
|
|
s = 1.0d0/dabs(z(k))
|
|
call dscal(n,s,z,1)
|
|
ynorm = s*ynorm
|
|
130 continue
|
|
140 continue
|
|
s = 1.0d0/dasum(n,z,1)
|
|
call dscal(n,s,z,1)
|
|
ynorm = s*ynorm
|
|
c
|
|
c solve u*z = w
|
|
c
|
|
do 160 kb = 1, n
|
|
k = n + 1 - kb
|
|
if (dabs(z(k)) .le. dabs(abd(m,k))) go to 150
|
|
s = dabs(abd(m,k))/dabs(z(k))
|
|
call dscal(n,s,z,1)
|
|
ynorm = s*ynorm
|
|
150 continue
|
|
if (abd(m,k) .ne. 0.0d0) z(k) = z(k)/abd(m,k)
|
|
if (abd(m,k) .eq. 0.0d0) z(k) = 1.0d0
|
|
lm = min0(k,m) - 1
|
|
la = m - lm
|
|
lz = k - lm
|
|
t = -z(k)
|
|
call daxpy(lm,t,abd(la,k),1,z(lz),1)
|
|
160 continue
|
|
c make znorm = 1.0
|
|
s = 1.0d0/dasum(n,z,1)
|
|
call dscal(n,s,z,1)
|
|
ynorm = s*ynorm
|
|
c
|
|
if (anorm .ne. 0.0d0) rcond = ynorm/anorm
|
|
if (anorm .eq. 0.0d0) rcond = 0.0d0
|
|
return
|
|
end
|
|
|
|
subroutine dgeco(a,lda,n,ipvt,rcond,z)
|
|
integer lda,n,ipvt(1)
|
|
double precision a(lda,1),z(1)
|
|
double precision rcond
|
|
c
|
|
c dgeco factors a double precision matrix by gaussian elimination
|
|
c and estimates the condition of the matrix.
|
|
c
|
|
c if rcond is not needed, dgefa is slightly faster.
|
|
c to solve a*x = b , follow dgeco by dgesl.
|
|
c to compute inverse(a)*c , follow dgeco by dgesl.
|
|
c to compute determinant(a) , follow dgeco by dgedi.
|
|
c to compute inverse(a) , follow dgeco by dgedi.
|
|
c
|
|
c on entry
|
|
c
|
|
c a double precision(lda, n)
|
|
c the matrix to be factored.
|
|
c
|
|
c lda integer
|
|
c the leading dimension of the array a .
|
|
c
|
|
c n integer
|
|
c the order of the matrix a .
|
|
c
|
|
c on return
|
|
c
|
|
c a an upper triangular matrix and the multipliers
|
|
c 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 rcond double precision
|
|
c an estimate of the reciprocal condition of a .
|
|
c for the system a*x = b , relative perturbations
|
|
c in a and b of size epsilon may cause
|
|
c relative perturbations in x of size epsilon/rcond .
|
|
c if rcond is so small that the logical expression
|
|
c 1.0 + rcond .eq. 1.0
|
|
c is true, then a may be singular to working
|
|
c precision. in particular, rcond is zero if
|
|
c exact singularity is detected or the estimate
|
|
c underflows.
|
|
c
|
|
c z double precision(n)
|
|
c a work vector whose contents are usually unimportant.
|
|
c if a is close to a singular matrix, then z is
|
|
c an approximate null vector in the sense that
|
|
c norm(a*z) = rcond*norm(a)*norm(z) .
|
|
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 linpack dgefa
|
|
c blas daxpy,ddot,dscal,dasum
|
|
c fortran dabs,dmax1,dsign
|
|
c
|
|
c internal variables
|
|
c
|
|
double precision ddot,ek,t,wk,wkm
|
|
double precision anorm,s,dasum,sm,ynorm
|
|
integer info,j,k,kb,kp1,l
|
|
c
|
|
c
|
|
c compute 1-norm of a
|
|
c
|
|
anorm = 0.0d0
|
|
do 10 j = 1, n
|
|
anorm = dmax1(anorm,dasum(n,a(1,j),1))
|
|
10 continue
|
|
c
|
|
c factor
|
|
c
|
|
call dgefa(a,lda,n,ipvt,info)
|
|
c
|
|
c rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) .
|
|
c estimate = norm(z)/norm(y) where a*z = y and trans(a)*y = e .
|
|
c trans(a) is the transpose of a . the components of e are
|
|
c chosen to cause maximum local growth in the elements of w where
|
|
c trans(u)*w = e . the vectors are frequently rescaled to avoid
|
|
c overflow.
|
|
c
|
|
c solve trans(u)*w = e
|
|
c
|
|
ek = 1.0d0
|
|
do 20 j = 1, n
|
|
z(j) = 0.0d0
|
|
20 continue
|
|
do 100 k = 1, n
|
|
if (z(k) .ne. 0.0d0) ek = dsign(ek,-z(k))
|
|
if (dabs(ek-z(k)) .le. dabs(a(k,k))) go to 30
|
|
s = dabs(a(k,k))/dabs(ek-z(k))
|
|
call dscal(n,s,z,1)
|
|
ek = s*ek
|
|
30 continue
|
|
wk = ek - z(k)
|
|
wkm = -ek - z(k)
|
|
s = dabs(wk)
|
|
sm = dabs(wkm)
|
|
if (a(k,k) .eq. 0.0d0) go to 40
|
|
wk = wk/a(k,k)
|
|
wkm = wkm/a(k,k)
|
|
go to 50
|
|
40 continue
|
|
wk = 1.0d0
|
|
wkm = 1.0d0
|
|
50 continue
|
|
kp1 = k + 1
|
|
if (kp1 .gt. n) go to 90
|
|
do 60 j = kp1, n
|
|
sm = sm + dabs(z(j)+wkm*a(k,j))
|
|
z(j) = z(j) + wk*a(k,j)
|
|
s = s + dabs(z(j))
|
|
60 continue
|
|
if (s .ge. sm) go to 80
|
|
t = wkm - wk
|
|
wk = wkm
|
|
do 70 j = kp1, n
|
|
z(j) = z(j) + t*a(k,j)
|
|
70 continue
|
|
80 continue
|
|
90 continue
|
|
z(k) = wk
|
|
100 continue
|
|
s = 1.0d0/dasum(n,z,1)
|
|
call dscal(n,s,z,1)
|
|
c
|
|
c solve trans(l)*y = w
|
|
c
|
|
do 120 kb = 1, n
|
|
k = n + 1 - kb
|
|
if (k .lt. n) z(k) = z(k) + ddot(n-k,a(k+1,k),1,z(k+1),1)
|
|
if (dabs(z(k)) .le. 1.0d0) go to 110
|
|
s = 1.0d0/dabs(z(k))
|
|
call dscal(n,s,z,1)
|
|
110 continue
|
|
l = ipvt(k)
|
|
t = z(l)
|
|
z(l) = z(k)
|
|
z(k) = t
|
|
120 continue
|
|
s = 1.0d0/dasum(n,z,1)
|
|
call dscal(n,s,z,1)
|
|
c
|
|
ynorm = 1.0d0
|
|
c
|
|
c solve l*v = y
|
|
c
|
|
do 140 k = 1, n
|
|
l = ipvt(k)
|
|
t = z(l)
|
|
z(l) = z(k)
|
|
z(k) = t
|
|
if (k .lt. n) call daxpy(n-k,t,a(k+1,k),1,z(k+1),1)
|
|
if (dabs(z(k)) .le. 1.0d0) go to 130
|
|
s = 1.0d0/dabs(z(k))
|
|
call dscal(n,s,z,1)
|
|
ynorm = s*ynorm
|
|
130 continue
|
|
140 continue
|
|
s = 1.0d0/dasum(n,z,1)
|
|
call dscal(n,s,z,1)
|
|
ynorm = s*ynorm
|
|
c
|
|
c solve u*z = v
|
|
c
|
|
do 160 kb = 1, n
|
|
k = n + 1 - kb
|
|
if (dabs(z(k)) .le. dabs(a(k,k))) go to 150
|
|
s = dabs(a(k,k))/dabs(z(k))
|
|
call dscal(n,s,z,1)
|
|
ynorm = s*ynorm
|
|
150 continue
|
|
if (a(k,k) .ne. 0.0d0) z(k) = z(k)/a(k,k)
|
|
if (a(k,k) .eq. 0.0d0) z(k) = 1.0d0
|
|
t = -z(k)
|
|
call daxpy(k-1,t,a(1,k),1,z(1),1)
|
|
160 continue
|
|
c make znorm = 1.0
|
|
s = 1.0d0/dasum(n,z,1)
|
|
call dscal(n,s,z,1)
|
|
ynorm = s*ynorm
|
|
c
|
|
if (anorm .ne. 0.0d0) rcond = ynorm/anorm
|
|
if (anorm .eq. 0.0d0) rcond = 0.0d0
|
|
return
|
|
end
|
|
|
|
|
|
subroutine dgedi(a,lda,n,ipvt,det,work,job)
|
|
integer lda,n,ipvt(1),job
|
|
double precision a(lda,1),det(2),work(1)
|
|
c
|
|
c dgedi computes the determinant and inverse of a matrix
|
|
c using the factors computed by dgeco or dgefa.
|
|
c
|
|
c on entry
|
|
c
|
|
c a double precision(lda, n)
|
|
c the output from dgeco or dgefa.
|
|
c
|
|
c lda integer
|
|
c the leading dimension of the array a .
|
|
c
|
|
c n integer
|
|
c the order of the matrix a .
|
|
c
|
|
c ipvt integer(n)
|
|
c the pivot vector from dgeco or dgefa.
|
|
c
|
|
c work double precision(n)
|
|
c work vector. contents destroyed.
|
|
c
|
|
c job integer
|
|
c = 11 both determinant and inverse.
|
|
c = 01 inverse only.
|
|
c = 10 determinant only.
|
|
c
|
|
c on return
|
|
c
|
|
c a inverse of original matrix if requested.
|
|
c otherwise unchanged.
|
|
c
|
|
c det double precision(2)
|
|
c determinant of original matrix if requested.
|
|
c otherwise not referenced.
|
|
c determinant = det(1) * 10.0**det(2)
|
|
c with 1.0 .le. dabs(det(1)) .lt. 10.0
|
|
c or det(1) .eq. 0.0 .
|
|
c
|
|
c error condition
|
|
c
|
|
c a division by zero will occur if the input factor contains
|
|
c a zero on the diagonal and the inverse is requested.
|
|
c it will not occur if the subroutines are called correctly
|
|
c and if dgeco has set rcond .gt. 0.0 or dgefa has set
|
|
c info .eq. 0 .
|
|
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,dswap
|
|
c fortran dabs,mod
|
|
c
|
|
c internal variables
|
|
c
|
|
double precision t
|
|
double precision ten
|
|
integer i,j,k,kb,kp1,l,nm1
|
|
c
|
|
c
|
|
c compute determinant
|
|
c
|
|
if (job/10 .eq. 0) go to 70
|
|
det(1) = 1.0d0
|
|
det(2) = 0.0d0
|
|
ten = 10.0d0
|
|
do 50 i = 1, n
|
|
if (ipvt(i) .ne. i) det(1) = -det(1)
|
|
det(1) = a(i,i)*det(1)
|
|
c ...exit
|
|
if (det(1) .eq. 0.0d0) go to 60
|
|
10 if (dabs(det(1)) .ge. 1.0d0) go to 20
|
|
det(1) = ten*det(1)
|
|
det(2) = det(2) - 1.0d0
|
|
go to 10
|
|
20 continue
|
|
30 if (dabs(det(1)) .lt. ten) go to 40
|
|
det(1) = det(1)/ten
|
|
det(2) = det(2) + 1.0d0
|
|
go to 30
|
|
40 continue
|
|
50 continue
|
|
60 continue
|
|
70 continue
|
|
c
|
|
c compute inverse(u)
|
|
c
|
|
if (mod(job,10) .eq. 0) go to 150
|
|
do 100 k = 1, n
|
|
a(k,k) = 1.0d0/a(k,k)
|
|
t = -a(k,k)
|
|
call dscal(k-1,t,a(1,k),1)
|
|
kp1 = k + 1
|
|
if (n .lt. kp1) go to 90
|
|
do 80 j = kp1, n
|
|
t = a(k,j)
|
|
a(k,j) = 0.0d0
|
|
call daxpy(k,t,a(1,k),1,a(1,j),1)
|
|
80 continue
|
|
90 continue
|
|
100 continue
|
|
c
|
|
c form inverse(u)*inverse(l)
|
|
c
|
|
nm1 = n - 1
|
|
if (nm1 .lt. 1) go to 140
|
|
do 130 kb = 1, nm1
|
|
k = n - kb
|
|
kp1 = k + 1
|
|
do 110 i = kp1, n
|
|
work(i) = a(i,k)
|
|
a(i,k) = 0.0d0
|
|
110 continue
|
|
do 120 j = kp1, n
|
|
t = work(j)
|
|
call daxpy(n,t,a(1,j),1,a(1,k),1)
|
|
120 continue
|
|
l = ipvt(k)
|
|
if (l .ne. k) call dswap(n,a(1,k),1,a(1,l),1)
|
|
130 continue
|
|
140 continue
|
|
150 continue
|
|
return
|
|
end
|
|
|
|
|