1      subroutine dgeco(a,lda,n,ipvt,rcond,z)
2      integer lda,n,ipvt(1)
3      double precision a(lda,1),z(1)
4      double precision rcond
5c
6c     dgeco factors a double precision matrix by gaussian elimination
7c     and estimates the condition of the matrix.
8c
9c     if  rcond  is not needed, dgefa is slightly faster.
10c     to solve  a*x = b , follow dgeco by dgesl.
11c     to compute  inverse(a)*c , follow dgeco by dgesl.
12c     to compute  determinant(a) , follow dgeco by dgedi.
13c     to compute  inverse(a) , follow dgeco by dgedi.
14c
15c     on entry
16c
17c        a       double precision(lda, n)
18c                the matrix to be factored.
19c
20c        lda     integer
21c                the leading dimension of the array  a .
22c
23c        n       integer
24c                the order of the matrix  a .
25c
26c     on return
27c
28c        a       an upper triangular matrix and the multipliers
29c                which were used to obtain it.
30c                the factorization can be written  a = l*u  where
31c                l  is a product of permutation and unit lower
32c                triangular matrices and  u  is upper triangular.
33c
34c        ipvt    integer(n)
35c                an integer vector of pivot indices.
36c
37c        rcond   double precision
38c                an estimate of the reciprocal condition of  a .
39c                for the system  a*x = b , relative perturbations
40c                in  a  and  b  of size  epsilon  may cause
41c                relative perturbations in  x  of size  epsilon/rcond .
42c                if  rcond  is so small that the logical expression
43c                           1.0 + rcond .eq. 1.0
44c                is true, then  a  may be singular to working
45c                precision.  in particular,  rcond  is zero  if
46c                exact singularity is detected or the estimate
47c                underflows.
48c
49c        z       double precision(n)
50c                a work vector whose contents are usually unimportant.
51c                if  a  is close to a singular matrix, then  z  is
52c                an approximate null vector in the sense that
53c                norm(a*z) = rcond*norm(a)*norm(z) .
54c
55c     linpack. this version dated 08/14/78 .
56c     cleve moler, university of new mexico, argonne national lab.
57c
58c     subroutines and functions
59c
60c     linpack dgefa
61c     blas daxpy,ddot,dscal,dasum
62c     fortran dabs,dmax1,dsign
63c
64c     internal variables
65c
66      double precision ddot,ek,t,wk,wkm
67      double precision anorm,s,dasum,sm,ynorm
68      integer info,j,k,kb,kp1,l
69c
70c
71c     compute 1-norm of a
72c
73      anorm = 0.0d0
74      do 10 j = 1, n
75         anorm = dmax1(anorm,dasum(n,a(1,j),1))
76   10 continue
77c
78c     factor
79c
80      call dgefa(a,lda,n,ipvt,info)
81c
82c     rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) .
83c     estimate = norm(z)/norm(y) where  a*z = y  and  trans(a)*y = e .
84c     trans(a)  is the transpose of a .  the components of  e  are
85c     chosen to cause maximum local growth in the elements of w  where
86c     trans(u)*w = e .  the vectors are frequently rescaled to avoid
87c     overflow.
88c
89c     solve trans(u)*w = e
90c
91      ek = 1.0d0
92      do 20 j = 1, n
93         z(j) = 0.0d0
94   20 continue
95      do 100 k = 1, n
96         if (z(k) .ne. 0.0d0) ek = dsign(ek,-z(k))
97         if (dabs(ek-z(k)) .le. dabs(a(k,k))) go to 30
98            s = dabs(a(k,k))/dabs(ek-z(k))
99            call dscal(n,s,z,1)
100            ek = s*ek
101   30    continue
102         wk = ek - z(k)
103         wkm = -ek - z(k)
104         s = dabs(wk)
105         sm = dabs(wkm)
106         if (a(k,k) .eq. 0.0d0) go to 40
107            wk = wk/a(k,k)
108            wkm = wkm/a(k,k)
109         go to 50
110   40    continue
111            wk = 1.0d0
112            wkm = 1.0d0
113   50    continue
114         kp1 = k + 1
115         if (kp1 .gt. n) go to 90
116            do 60 j = kp1, n
117               sm = sm + dabs(z(j)+wkm*a(k,j))
118               z(j) = z(j) + wk*a(k,j)
119               s = s + dabs(z(j))
120   60       continue
121            if (s .ge. sm) go to 80
122               t = wkm - wk
123               wk = wkm
124               do 70 j = kp1, n
125                  z(j) = z(j) + t*a(k,j)
126   70          continue
127   80       continue
128   90    continue
129         z(k) = wk
130  100 continue
131      s = 1.0d0/dasum(n,z,1)
132      call dscal(n,s,z,1)
133c
134c     solve trans(l)*y = w
135c
136      do 120 kb = 1, n
137         k = n + 1 - kb
138         if (k .lt. n) z(k) = z(k) + ddot(n-k,a(k+1,k),1,z(k+1),1)
139         if (dabs(z(k)) .le. 1.0d0) go to 110
140            s = 1.0d0/dabs(z(k))
141            call dscal(n,s,z,1)
142  110    continue
143         l = ipvt(k)
144         t = z(l)
145         z(l) = z(k)
146         z(k) = t
147  120 continue
148      s = 1.0d0/dasum(n,z,1)
149      call dscal(n,s,z,1)
150c
151      ynorm = 1.0d0
152c
153c     solve l*v = y
154c
155      do 140 k = 1, n
156         l = ipvt(k)
157         t = z(l)
158         z(l) = z(k)
159         z(k) = t
160         if (k .lt. n) call daxpy(n-k,t,a(k+1,k),1,z(k+1),1)
161         if (dabs(z(k)) .le. 1.0d0) go to 130
162            s = 1.0d0/dabs(z(k))
163            call dscal(n,s,z,1)
164            ynorm = s*ynorm
165  130    continue
166  140 continue
167      s = 1.0d0/dasum(n,z,1)
168      call dscal(n,s,z,1)
169      ynorm = s*ynorm
170c
171c     solve  u*z = v
172c
173      do 160 kb = 1, n
174         k = n + 1 - kb
175         if (dabs(z(k)) .le. dabs(a(k,k))) go to 150
176            s = dabs(a(k,k))/dabs(z(k))
177            call dscal(n,s,z,1)
178            ynorm = s*ynorm
179  150    continue
180         if (a(k,k) .ne. 0.0d0) z(k) = z(k)/a(k,k)
181         if (a(k,k) .eq. 0.0d0) z(k) = 1.0d0
182         t = -z(k)
183         call daxpy(k-1,t,a(1,k),1,z(1),1)
184  160 continue
185c     make znorm = 1.0
186      s = 1.0d0/dasum(n,z,1)
187      call dscal(n,s,z,1)
188      ynorm = s*ynorm
189c
190      if (anorm .ne. 0.0d0) rcond = ynorm/anorm
191      if (anorm .eq. 0.0d0) rcond = 0.0d0
192      return
193      end
194