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