subroutine decomp(ndim,n,a,cond,ipvt,work)
c
integer ndim,n
double precision a(ndim,n),cond,work(n)
integer ipvt(n)
c
c decomposes a double precision matrix by gaussian elimination
c and estimates the condition of the matrix.
c
c use solve to compute solutions to linear systems.
c
c input..
c
c ndim = declared row dimension of the array containing a.
c
c n = order of the matrix.
c
c a = matrix to be triangularized.
c
c output..
c
c a contains an upper triangular matrix u and a permuted
c version of a lower triangular matrix i-l so that
c (permutation matrix)*a = l*u .
c
c cond = an estimate of the condition of a .
c for the linear system a*x = b, changes in a and b
c may cause changes cond times as large in x .
c if cond+1.0 .eq. cond , a is singular to working
c precision. cond is set to 1.0d+32 if exact
c singularity is detected.
c
c ipvt = the pivot vector.
c ipvt(k) = the index of the k-th pivot row
c ipvt(n) = (-1)**(number of interchanges)
c
c work space.. the vector work must be declared and included
c in the call. its input contents are ignored.
c its output contents are usually unimportant.
c
c the determinant of a can be obtained on output by
c det(a) = ipvt(n) * a(1,1) * a(2,2) * ... * a(n,n).
c
double precision ek, t, anorm, ynorm, znorm
integer nm1, i, j, k, kp1, kb, km1, m
double precision dabs, dsign
c
ipvt(n) = 1
if (n .eq. 1) go to 80
nm1 = n - 1
c
c compute 1-norm of a
c
anorm = 0.0d0
do 10 j = 1, n
t = 0.0d0
do 5 i = 1, n
t = t + dabs(a(i,j))
5 continue
if (t .gt. anorm) anorm = t
10 continue
c
c gaussian elimination with partial pivoting
c
do 35 k = 1,nm1
kp1= k+1
c
c find pivot
c
m = k
do 15 i = kp1,n
if (dabs(a(i,k)) .gt. dabs(a(m,k))) m = i
15 continue
ipvt(k) = m
if (m .ne. k) ipvt(n) = -ipvt(n)
t = a(m,k)
a(m,k) = a(k,k)
a(k,k) = t
c
c skip step if pivot is zero
c
if (t .eq. 0.0d0) go to 35
c
c compute multipliers
c
do 20 i = kp1,n
a(i,k) = -a(i,k)/t
20 continue
c
c interchange and eliminate by columns
c
do 30 j = kp1,n
t = a(m,j)
a(m,j) = a(k,j)
a(k,j) = t
if (t .eq. 0.0d0) go to 30
do 25 i = kp1,n
a(i,j) = a(i,j) + a(i,k)*t
25 continue
30 continue
35 continue
c
c cond = (1-norm of a)*(an estimate of 1-norm of a-inverse)
c estimate obtained by one step of inverse iteration for the
c small singular vector. this involves solving two systems
c of equations, (a-transpose)*y = e and a*z = y where e
c is a vector of +1 or -1 chosen to cause growth in y.
c estimate = (1-norm of z)/(1-norm of y)
c
c solve (a-transpose)*y = e
c
do 50 k = 1, n
t = 0.0d0
if (k .eq. 1) go to 45
km1 = k-1
do 40 i = 1, km1
t = t + a(i,k)*work(i)
40 continue
45 ek = 1.0d0
if (t .lt. 0.0d0) ek = -1.0d0
if (a(k,k) .eq. 0.0d0) go to 90
work(k) = -(ek + t)/a(k,k)
50 continue
do 60 kb = 1, nm1
k = n - kb
t = 0.0d0
kp1 = k+1
do 55 i = kp1, n
t = t + a(i,k)*work(i)
55 continue
work(k) = t + work(k)
m = ipvt(k)
if (m .eq. k) go to 60
t = work(m)
work(m) = work(k)
work(k) = t
60 continue
c
ynorm = 0.0d0
do 65 i = 1, n
ynorm = ynorm + dabs(work(i))
65 continue
c
c solve a*z = y
c
call solve(ndim, n, a, work, ipvt)
c
znorm = 0.0d0
do 70 i = 1, n
znorm = znorm + dabs(work(i))
70 continue
c
c estimate condition
c
cond = anorm*znorm/ynorm
if (cond .lt. 1.0d0) cond = 1.0d0
return
c
c 1-by-1
c
80 cond = 1.0d0
if (a(1,1) .ne. 0.0d0) return
c
c exact singularity
c
90 cond = 1.0d+32
return
end